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Abstract 



In this article we study both left-invariant (convection-)diffusions and left-invariant Hamilton- Jacobi equations 
(erosions) on the space x of 3D-positions and orientations naturally embedded in the group SE{3) of 3D- 
rigid body movements. The general motivation for these (convection-)diffusions and erosions is to obtain crossing- 
preserving fiber enhancement on probability densities defined on the space of positions and orientations. 

The linear left-invariant (convection-)diffusions are forward Kolmogorov equations of Brownian motions on xi 
S'^ and can be solved by x 5*^ -convolution with the corresponding Green's functions or by a finite difference 
scheme. The left-invariant Hamilton- Jacobi equations are Bellman equations of cost processes on R^ x S'^ and they 
are solved by a morphological R^ x ^S*^ -convolution with the corresponding Green's functions. We will reveal the 
remarkable analogy between these erosions/dilations and diffusions. Furthermore, we consider pseudo-linear scale 
spaces on the space of positions and orientations that combines dilation and diffusion in a single evolution. 

In our design and analysis for appropriate linear, non-linear, morphological and pseudo-linear scale spaces on 
R^ X S'^ we employ the underlying differential geometry on SE{3), where the frame of left-invariant vector fields 
serves as a moving frame of reference. Furthermore, we will present new and simpler finite difference schemes for 
our diffusions, which are clear improvements of our previous finite difference schemes. 

We apply our theory to the enhancement of fibres in magnetic resonance imaging (MRI) techniques (HARDI 
and DTI) for imaging water diffusion processes in fibrous tissues such as brain white matter and muscles. We pro- 
vide experiments of our crossing-preserving (non-linear) left-invariant evolutions on neural images of a human brain 
containing crossing fibers. 

Keywords: Nonlinear diffusion, Lie groups, Hamilton- Jacobi equations, Partial differential equations, Sub-Riemannian geometry, 
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1 Introduction 



Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI) involves magnetic resonance techniques for non-invasively 
measuring local water diffusion inside tissue. Local water diffusion profiles reflect underlying biological fiber structure 
of the imaged area. For instance in the brain, diffusion is less constrained parallel to nerve fibers than perpendicular 
to them. In this way, the measurement of water diffusion gives information about the fiber structures present, which 
allows extraction of clinical information from these scans. 



The diffusion of water molecules in tissue over time is described by a transition density y ^ PtiYt = J \ = Jo) 
that reflects the probability density of finding a water molecule at time t > at position y G given that it started 
at y G at time t = 0. Here the family of random variables {Yt)t>o (stochastic process) with joint state space 
reflects the distribution of water molecules over time. The function pt can be directly related to MRI signal attenuation 
of Diffusion- Weighted image sequences, so can be estimated given enough measurements. The exact methods to do 
this are described by e.g. Alexander |3|. Diffusion tensor imaging (DTI), introduced by Basser et al. ifTTIl . assumes 
that Pt can be described for each position y G by an anisotropic Gaussian function. So 

/ I X 1 (y^-y)^(I^(y))-l(y^-y) 

Pt{Yt = y \Yo=y) = ----^^===e 
(47rt)2^|detL)(y)| 

where is a tensor field of positive definite symmetric tensors on R^. In a DTI- visualization one plots the surfaces 

y + {vGR^ I v^I)-i(y)v = /i2}, (1) 

where /i > is fixed and y e ft with ft some compact subset of R^. From now on we refer to these ellipsoidal surfaces 
as DTI-glyphs. The drawback of this anisotropic Gaussian approximation in DTI is the limited angular resolution of 
the corresponding probability density U y\ S'^ ^ R+ on positions and orientations 

and thereby unprocessed DTI is not capable of representing crossing, kissing or diverging fibers, cf. O. 



High Angular Resolution Diffusion Imaging (HARDI) is another recent magnetic resonance imaging technique for 
imaging water diffusion processes in fibrous tissues such as brain white matter and muscles. HARDI provides for each 
position in R^ and for each orientation an MRI signal attenuation profile, which can be related to the local diffusivity of 
water molecules in the corresponding direction. Typically, in HARDI modeling the Fourier transform of the estimated 
transition densities is considered at a fixed characteristic radius (generally known as the b-value), |23|. As a result, 
HARDI images are distributions (y, n) U (y, n) over positions and orientations, which are often normalized per 
position. HARDI is not restricted to functions on the 2-sphere induced by a quadratic form and is capable of reflecting 
crossing information, see Figure [T] where we visualize HARDI by glyph visualization as defined below. 



Definition 1 A glyph of a distribution U : R^ x S*^ ^ R+ on positions and orientations is a surface S^{U)(y) = 
{y -\- fiU{y^n) n \ n G S'^} C for some j G R^ and /i > 0. A glyph visualization of the distribution 
U : R^ X S*^ ^ R+ is a visualization of a field y ^ S^{U){y) of glyphs, where fi > is a suitable constant. 



For the purpose of detection and visualization of biological fibers, DTI and HARDI data enhancement should maintain 
fiber junctions and crossings, while reducing high frequency noise in the joined domain of positions and orientations. 
Promising research has been done on constructing diffusion/regularization processes on the 2-sphere defined at each 
spatial locus separately ||22l[37l[38l[66l as an essential pre-processing step for robust fiber tracking. In these approaches 
position and orientation space are decoupled, and diffusion is only performed over the angular part, disregarding spatial 
context. Consequently, these methods are inadequate for spatial denoising and enhancement, and tend to fail precisely 
at the interesting locations where fibres cross or bifurcate. 
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Figure 1 : This figure shows glyph visuaHzations of HARDI and DTI-images of a 2D-sHce in the brain where neural 
fibers in the corona radiata cross with neural fibers in the corpus callosum. Here DTI and HARDI are visualized 
differently; HARDI is visualized according to Def.[T] whereas DTI is visualized using Eq. ([T]). 



In contrast to previous work on diffusion of DW-MRI |l22l |37l |38l [66l [Sjl, we consider both the spatial and the 
orientational part to be included in the domain, so a HARDI dataset is considered as a function U : x S*^ ^ R. 
Furthermore, we explicitly employ the proper underlying group structure, that naturally arises by embedding the 
coupled space of positions and orientations 

R^ y\ := SE{3)/{{0} x 50(2)), 

as the quotient of left cosets, into the group SE{3) = R^ x 5*0 (3) of 3D-rigid motions. The relevance of group theory 
in DTI/HARDI (DW-MRI) imaging has also been stressed in promising and well-founded recent works 14411451 1461 . 
However these works rely on bi-invariant Riemannian metrics on compact groups (such as 6*0 (3)) and in our case the 
group SE{3) is neither compact nor does it permit a bi-invariant metric |[6ll3T1l. 

Throughout this article we use the following identification between the DW-MRI data (y, n) /7(y, n) and functions 
U : 5'^(3) ^R given by 

U{y^R) = U{y,Re,). (3) 

By definition one has U{y^ RR^^^a) = U{y^R) for all a G [0,27r), where Re^,a is the counterclockwise rotation 
around by angle a. 

In general the advantage of our approach on ^£^(3) is that we can enhance the original HARDI/DTI data using 
simultaneously orientational and spatial neighborhood information, which potentially leads to improved enhancement 
and detection algorithms, | 34, 60, 59 1. See Figure [2] where fast practical implementations |60 | of the theory developed 
in [ 34, ch:8.2] have been applied. The potential clinical impact is evident: By hypo-elliptic diffusions on R^ x 5*^ one 
can generate distributions from DTI that are similar to HARDI-modeling as recently reported by Prckovska et al. (591 . 
This allows a reduction of scanning directions in areas where the random walks processes that underly hypo-elliptic 
diffusion [34, ch:4.2] on R^ x 5*^ yield reasonable fiber extrapolations. Experiments on neural DW-MRI images 
containing crossing fibers of the corpus callosum and corona radiata show that extrapolation of DTI (via hypo-elliptic 
diffusion) can cope with HARDI for a whole range of reasonable 6- values, |59|. See Figure [2] However, on the 
locations of crossings HARDI in principle produces more detailed information than extrapolated DTI and application 
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Figure 2: DTI and HARDI data containing fibers of the corpus callosum and the corona radiata in a human brain with b- 
value lOOOs/mm^ on voxels of (2mm)^, cf. |59|. We visuaHze a 10 x 16-sHce (162 samples on S'^ using icosahedron 
tessellations) of interest from 104 x 104 x 10 x (162 x 3) datasets. Top row, region of interest with fractional anisotropy 
intensities with colorcoded DTI-principal directions. Middle row, DTI data U (49 scanned orientations) vizualized 
according to Eq. ([T]) respectively Def. [T] Bottom row: HARDI data (Q-ball with / < 4, (231) of the same region, 
processed DTI data ^t{U). For the sake of visualization we applied min-max normalization of n 1-^ $^([/)(y, n) for 
all positions y. For parameter settings of the hypo-elliptic diffusion operator see Section [6j Eq. ( [49| ). 



of the same hypo-elliptic diffusion on HARDI removes spurious crossings that arise in HARDI, see Figure |3] and the 
recent work |60|. 

In this article we will build on the recent previous work |[59j[34j[60l, and we address the following open issues that 
arise immediately from these works: 



• Can we adapt the diffusion on x S'^ locally to the initial HARDI/DTI image? 

• Can we apply left-invariant Hamilton- Jacobi equations (erosions) to sharpen the data without grey-scale trans- 
formations (squaring) needed in our previous work ? 

• Can we classify the viscosity solutions of these left-invariant Hamilton- Jacobi equations? 

• Can we find analytic approximations for viscosity solutions of these left-invariant Hamilton- Jacobi equations on 

X 5*^, likewise the analytic approximations we derived for linear left-invariant diffusions, cf. | [34| ch:6.2]? 

• Can we relate alternations of greyscale transformations and linear diffusions to alternations of linear diffusions 
and erosions? 
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Figure 3: Same settings as Fig: 2, except for different 6- value and a different region of interest and the hypo-elliptic 



diffusion, Eq. (49), is applied to the HARDI dataset. 



• The resolvent Green's functions of the direction process on y\ S'^ contain singularities at the origin, (341 
ch: 6. 1. l,Fig. 8]. Can we overcome this complication in our algorithms for the direction process and can we 
analyze iterations of multiple time integrated direction process to control the filling of gaps? 

• Can we combine left-invariant diffusions and left-invariant dilations/erosions in a single pseudo-linear scale 
space on x 5*^, generalizing the approach by Florack et al. | 36| for greyscale images to DW-MRI images? 



• Can we avoid spherical harmonic transforms and the sensitive regularization parameter treg IMl ch:7. 1,7.2] in 
our finite difference schemes and obtain both faster and simpler numerical approximations of the left-invariant 
vector fields ? 



To address these issues, we introduce besides linear scale spaces, morphological scale spaces and pseudo-linear scale 
spaces (y, n, t) ^ W(y, n, t), for all y G R^, n G S'^.t > 0, defined on (R^ x S'^) x R+, where we use the input 
DW-MRI image as initial condition iy(y, n, 0) = /7 (y, n). 

To get a preview of how these approaches perform on the same neural DTI dataset (different slice) considered in 1591 , 
see Fig. |5] where we used 

V(r/)(y,n) = f ,^^^'"\~^r"?v\ )'' = inWC^(y,n) I n £ S'}. (4) 



Typically, if linear diffusions are directly applied to DTI the fibers visible in DTI are propagated in too many directions. 
Therefore we combined these diffusions with monotonic transformations in the codomain R+, such as squaring input 
and output cf. 1 34 , ,60, ,59] . Visually, this produces anatomically plausible results, cf. Fig.[2]and Fig.|3j but does not 
allow large global variations in the data. This is often problematic around ventricle areas in the brain, where the 
glyphs are usually larger than those along the fibers, as can be seen in the top row of Fig. [5] In order to achieve a 
better way of sharpening the data where global maxima do not dominate the sharpening of the data, cf. Fig.[4j we 
propose morphological scale spaces on R^ x S'^ where transport takes place orthogonal to the fibers, both spatially 
and spherically, see Fig.|6] The result of such an erosion after application of a linear diffusion is depicted down left in 
Fig.|5] where the diffusion has created crossings in the fibers and where the erosion has visually sharpened the fibers. 
Simultaneous dilation and diffusion can be achieved in a pseudo-linear scale space that conjugates a diffusion with a 
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Figure 4: From left to right. Noisy artificial dataset, output diffused dataset (thresholded), squared output diffused 
dataset as in (591 [601, 5^ -eroded output, Eq.([52]), diffused dataset, Eq.(|49). 



specific grey-value transformation. An experiment of applying such a left-invariant pseudo-linear scale space to DTI- 
data is given up-right in Figure|5] Regarding the numerics of the evolutions we mainly consider the left-invariant finite 
difference approach in [ 34, ch:7], as an alternative to the analytic kernel implementations in [ 34, ch:8.2] and l59l[6Qll. 
Here we avoid the discrete spherical Harmonic transforms 1341 . but use fast precomputed linear interpolations instead 
as our finite difference schemes are of first order accuracy anyway. Regarding the Hamilton- Jacobi equations involved 
in morphological and pseudo-linear scale spaces we have, akin to the linear left-invariant diffusions | 34, ch:7], two 
options: analytic morphological x 5^ -convolutions and finite differences. Regarding fast computation on sparse 
grids the second approach is preferable. Regarding geometric analysis the first approach is preferable. 



We show that our morphological R^ x 5^ -convolutions with analytical morphological Green's functions are the unique 
viscosity solutions of the corresponding Hamilton- Jacobi equations on R^ x S*^ . Thereby, we generalize the results 



in [ 35, ch:10], [53 J (on the Heisenberg group) to Hamilton- Jacobi equations on the space 
orientations. 



S of positions and 



Evolutions on HARDI-DTI must commute with all rotations and translations. Therefore evolutions on HARDI and 
DTI and underlying metric (tensors) are expressed in a local frame of reference attached to fiber fragments. This frame 
{^1, . . . , ^e} consists of 6 left-invariant vector fields on SE{3), given by 



AiU{y,R) = lim 



U{{y,R)e^^^)-U{{y,R)e-^^^) 
2h 



(5) 



where {Ai, . . . , Ae} is the basis for the Lie-algebra at the unity element and Te{SE{3)) 3 A e SE{3) is the 

exponential map in SE{3) and where the group product on SE{3) is given by 

(x, R) (x^ R') = (x + Rx\ RR'), (6) 

for all positions x,x' G R^ and rotations R,R' G SO (3). The details will follow in Section |4j see also (341 
ch:3.3,Eq. 23-25] and t34, ch:7] for implementation. In order to provide a relevant intuitive preview of this moving 
frame of reference we refer to Fig. [6] The associated left-invariant dual frame {d^^, . . . , d^^} is uniquely determined 
by 

{dA\ Aj) := dA'iAj) = 5], i,i = 1, . . . , 6, (7) 
where (5' = 1 if i = j and zero else. Then all possible left-invariant metrics are given by 



G 
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d^* L ^ , (g) d^-^ L ^ , . 



(8) 



where y G R^, n e S'^ and where Rn G S0{3) is any rotation that maps = (0, 0, 1)^ onto the normal n G 5^, i.e. 



Rn^z n. 



(9) 



and where gij G C. Necessary and sufficient conditions on gij to induce a well-defined left-invariant metric on the 
quotient R^ >^ 6'^ = {SE{3)/{{0} x S0{2))) can be found in Appendix|E] It turns out that the matrix [gij] must be 
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Figure 5: DTI data of corpus callosum and corona radiata fibers in a human brain with 6- value 1 000s /mm^ on voxels 
of (2mm)^. Top row: DTI-visuaHzation according to Eq. ([l]). The yellow box contains 13 x 22 x 10 glyphs with 
162 orientations of the input DTI-data depicted in the left image of the middle row. This input-DTI image U is 
visualized using Eq. ^ and Rician noise rjr 1 34, Eq. 90] with a = 10 ^ has been included. Operator V is defined in 
Eq. (|4]). Middle row, right: output of finite diffe rence scheme pseudo-linear scale space (for parameters see Section 
\\a\ . Bottom row, left: output erosion, Eq. ( 136 ) after hypo-elliptic diffusion, Eq. (49 ) with {D^^ = 0.04, D^^ = 1, 
t = 1), right: output of non-linear adaptive diffusions discussed in [21 1 with Z)^^ = 0.01. All evolutions commute 
with rotations and translations and are implemented by finite differences (Section [12]) with time t and stepsize At. 
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Figure 6: A curve [0, 1] 3 s 7(5) = (x(s),n(s)) ^ >^ consists of a spatial part s x(s) (left) and an angular 
part s 1-^ n(s) (right). Along this curve we have the moving frame of reference {Ai\~(^^^^}i=l with 7(5) = (x(s), i^n(s)) where 
^n(s) ^ SO(3) is any rotation such that Rji(s)^z = n(s) G S^. Here Ai, with = A|(o i) denote the left-invariant vector fields 
in SE{3). To ensure that the diffusions and erosions do not depend on the choice Rn(s) ^ SO (3), Eq. |9j, these left-invariant 
evolution equations must be isotropic in the tangent planes span{^i , A2} and span{^4, ^5}. Diffusion/convection primarily takes 
place along A3 in space and (outward) in the plane span{^4, ^5} tangent to S^. Erosion takes place both inward in the tangent 
plane span{^i, ^2} in space and inward in the plane span{^4, ^s}- 



constant and diagonal gij = ^Sij, z, j = 1 . . . , 6 with D'' G M+ U 00, with D^^ = L)^^ D^^ = D^^ , D^^ = 0. 
Consequently, the metric is parameterized by the values D^^ ^ D^^ ^ D^^ and in the sequel we write 

G^''^''^"" := -^{dA^ ^ dA^ + dA^ dA^) + 7^(d^^ d^^) + -^{dA"" ^ dA"" + d^^ d^^) 

JJLl JJ66 JJ^^ 

The corresponding metric tensor on the quotient x S*^ = {SE{3) / ({0} x 50(2))) is given by 

Cr^^"(i M(y^r.),id^ A|(,.„)) = GgW(E ^.|(,^„) A-|(,,„)) (10) 

It is well-defined on x S"^ as the choice of i?n, Eq J9](, does not matter (right-multiplication with i?e^,a boils down 
to rotations in the isotropic planes depicted in Figure |6]i and with the differential operators on 9? x S"^: 

( A-I(,,n) U){y.n) = lim ^(y+^^ne..n)-^^(y-.i.„e..,n) ^ 

(^3+j|(y,n) U){y,n) = lim ^ ^ ^ , j = 1,2,3, 

where Re^^h denotes the counter-clockwise rotation around axis ej by angle h, with ei = (1,0,0), 62 = (0,1,0), 
63 = (0, 0, 1), which (except for ^3) do depend on the choice of R^ satisfying Eq. ([9]). 



In |[34j ch:6.2] we have analytically approximated the hypo-elliptic diffusion kernels for both the direction process and 
Brownian motion on the sub-Riemannian manifold (or contact manifold 1 13 1) {SE{3) , d^^ , dA'^ , d^^) using contrac- 
tion towards a nilpotent group. For the erosions we employ a similar type of technique to analytically approximate the 
erosion (and dilation) kernels that describe the growth of balls in the sub-Riemannian manifold {SE{3)^ dA^^ dA^). 

A sub-Riemannian manifold is a Riemannian manifold with the extra constraint that tangent-vectors to curves in 
that Riemannian manifold are not allowed to use certain subspaces of the tangent space. For example, curves in 
{SE{3),dA^,dA'^, dA^) are curves 7 : [0, L] SE{3) with the constraint 

{dA\^^A{s)) = {dA\^^A{s)) = {dA%^^^ A{s)) = 0, (11) 

for all s e [0, L], L > 0. Note that Eq. ([TT]) implies that 

^(s) = {dA\^^ Ais)) ^3|^(,) + {dA\^^A{s)) A^\~^^^ + {dA'\^^^^ A{s)) A,\~^^^ 
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Curves satisfying ( pTj ) are called horizontal curves and we visualized a horizontal curve in Figure [6] For further details 
on differential geometry in (sub-Riemannian manifolds within) SE(?t) and x S'^ see Appendices [a| [b| and [c{ [e| 
andEl 



1.1 Outline of the article 



This paper is organized as follows. In section [2| w e will explain the embedding of the coupled space of positions 
and orientations x S'^ into SE{3). In sectionplwe explain why operators on DW-MRI must be left-invariant and 
we consider, as an example, convolutions on x S'^. In section |4 we construct the left-invariant vector fields on 
SE{3). In general there are two ways of considering vector fields. Either one considers them as differential operators 
on smooth locally defined functions, or one considers them as tangent vectors to equivalent classes of curves. These 
two viewpoints are equivalent, for a formal proof see | 7, Prop. 2.4]. Throughout this article we will consider them as 
differential operators and use them as reference frame, cf. Fig|6] in our evolution equations in later sections. Then in 
Sectionjsjwe consider morphological convolutions on R^ x 5^. 

In Section [6] we consider all possible linear left-invariant diffusions that are solved by R^ x 5^ -convolution (Section 
[3| with the corresponding Green's function. Subsequently, in Section|7]we consider their morphological counter part: 
left-invariant Hamilton- Jacobi equations (i.e. erosions), the viscosity solutions of which are given by morphologi- 
cal convolution, cf. Section [5] with the corresponding (morphological) Green's function. This latter result is a new 
fundamental mathematical result, a detailed proof is given in Appendix |B] 

Subsequently, in Section [8] and in Section [9] we provide respectively the underlying probability theory and the un- 
derlying Cartan differential geometry of the evolutions considered in Section [5] and Section [6] Then in Section [10 
we derive analytic approximations for the Green's functions of both the linear and the morphological evolutions on 
R^ X 6*^. Section [IT] deals with pseudo linear scale spaces which are evolutions that combine the non-linear generator 
of erosions with the generator of diffusion in a single generator. 



Section 12 deals with the numerics of the (convection)-diffusions, the Hamilton- Jacobi equations and the pseudo linear 



scale spaces evolutions. Section 13 summarizes our work on adaptive left-invariant diffusions on DW-MRI data, for 



more details we refer to the Master thesis by Eric Creusen 121 J . Finally, we provide experiments in Section 14 



2 The Embedding of >^ 5^2 ^^^^ 5<^(^3) 



In order to generalize our previous work on line/contour-enhancement via left-invariant diffusions on invertible orien- 
tation scores of 2I)-images we first investigate the group structure on the domain of an HARDI image. Just like orien- 
tation scores are scalar-valued functions on the coupled space of 2D-positions and orientations, i.e. the 2I)-Euclidean 
motion group, HARDI images are scalar- valued functions on the coupled space of 3D-positions and orientations. This 
generalization involves some technicalities since the 2-sphere S'^ = {x G R^ | ||x|| = 1} is not a Lie-group prope][^ 
in contrast to the 1-sphere — {x G R^ | ||x|| = 1}. To overcome this problem we will embed R^ x into 
5'^(3) which is the group of 3I)-rotations and translations (i.e. the group of 3D-rigid motions). As a concatenation 
of two rigid body-movements is again a rigid body movement, the product on SE(Z) is given by ([6]). The group 
SE(Z) is a semi-direct product of the translation group R^ and the rotation group 5*0 (3), since it uses an isomorphism 
1-^ (x i?x) from the rotation group onto the automorphisms on R^. Therefore we write R^ x 5*0 (3) rather than 
R^ X S'0(3) which would yield a direct product. The groups SEiZ) and 6'0(3) are not commutative. Throughout 
this article we will use Euler- angle parametrization for 5*0 (3), i.e. we write every rotation as a product of a rotation 

^If were a Lie-group then its left-invariant vector fields would be non-zero everywhere, contradicting Poincare's "hairy ball theorem" (proven 
by Brouwer in 1912), or more generally the Poincare-Hopf theorem (the Euler-characteristic of an even dimensional sphere S'^^ is 2). 
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around the 2:-axis, a rotation around the ?/-axis and a rotation around the z-axis again. 

R = Re^,-fRey,(3Re^,a , (12) 

where all rotations are counter-clockwise. ExpHcit formulas for matrices R^^ , R^^ , R^^ ^q,, are given in (391 ch:7.3. 1]. 
The advantage of the Euler angle parametrization is that it directly parameterizes S0{^) / S0{2) = 5*^ as well. 
Here we recall that S0{^) / S0{2) denotes the partition of all left cosets which are equivalence classes [g] = {h e 
SO{3) \ h ^ g} = g S0{2) under the equivalence relation gi ^ g2 <^ 9i^92 ^ S0{2) where we identified S0{2) 
with rotations around the z-axis and we have 

SO{3)/SO{2) 9 [i?e,,7^e„/3] = {i?e.,7^e„,/3i?e.,a \ a G [0, 27r)} O 

n(/3,7) := (cos 7 sin sin 7 sin /3, cos (3) = Re^^^Rey^(3Re^,a^z ^ S'^. 

Like all parameterizations of SO{3)/SO{2), the Euler angle parametrization suffers from the problem that there does 
not exist a global diffeomorphism from a sphere to a plane. In the Euler-angle parametrization the ambiguity arises at 
the north and south-poles: 

Rez,-r Rey,p^O Rez,oi = i^e^ ,7-5^6^ ,/3=0^e^ ,a+(5 , and i?e^ ,7 ,/3=7r ,a = Re^ R^y R^z ,ci-\-S , for all ^ G [0, 27r) . 

(14) 

Consequently, we occasionally need a second chart to cover 50(3); 

R = Re^,^ReyjRe^,a, (15) 

which again implicitly parameterizes 5*0 (3) /5'0 (2) = 5^ using different ball-coordinates /5 G [— 7r,7r),7 G (— |, |), 
n(;^,7) = i?e,,7^e^,^e^ (sin^, -cos;^ sin7,coS;^ cos7)^, (16) 
but which has ambiguities at the intersection of the equator with the x-axis, 1341 . 

Re,,7\,^=±f^^z,a = ^e,,7T5^e^,^=±f ^e,,a±5, for all (5 G [0, 27r) , (17) 

see Figure [T] Away from the intersection of the z and x-axis with the sphere one can accomplish conversion between 
the two charts by solving for for either (a, ^, 7) or (a, /3, 7) in R^^^^R^^^R^^^a = i^e,,7^e^,/3^e,,a. 

Now that we have explained the isomorphism n = Rez G 5^ ^ SO{3)/SO{2) 3 [R] explicitly in charts, we return 
to the domain of HARDI images. Considered as a set, this domain equals the space of 3D-positions and orientations 
X S*^. However, in order to stress the fundamental embedding of the HARDI-domain in ^£^(3) and the thereby 
induced (quotient) group- structure we write x S*^, which is given by the following Lie-group quotient: 

R^ X := (R^ X 5O(3))/({0} x 50(2)). 

Here the equivalence relation on the group of rigid-motions ^£^(3) = R^ x 50(3) equals 

(x, R) - (x', R') o (x, R)~^{x\ R') G {0} X 50(2) ^ x = x' and R'^R' is a rotation around z-axis 

and set of equivalence classes within 5£^(3) under this equivalence relation (i.e. left cosets) equals the space of 
coupled orientations and positions and is denoted by R^ x 5^ . 



3 Linear Convolutions on xi S'^ 



In this article we will consider convection-diffusion operators on the space of HARDI images. We shall model the 
space of HARDI images by the space of square integrable functions on the coupled space of positions and orientations, 
i.e. L2 (R^ X 5^). We will first show that such operators should be left-invariant with respect to the left-action of SE{3) 
onto the space of HARDI images. This left-action of the group 5£^(3) onto R^ x 5^ is given by 

g ■ (y, n) = {Ry + x, i?n), g = (x, R) G 5^(3), x, y G R^ n G 5^ G 50(3) (18) 

and it induces the so-called left-regular action of the same group on the space of HARDI images similar to the left- 
regular action on 3D-images (for example orientation-marginals of HARDI images): 
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Figure 7: The two charts which together appropriately parameterize the sphere 5*^ = SO(Z) j S0(2) where the 



rotation-parameters ot and a are free. The first chart (left-image) is the common Euler- angle parametrization ( 12 ), the 



second chart is given by ( 15 ). The first chart has singularities at the north and at the south-pole (inducing ill-defined 



parametrization of the left-invariant vector fields ( 3 1 ) at the unity element) whereas the second chart has singularities 
at (±1,0,0). 



Definition 2 The left-regular actions ofSE{3) onto L2(M^ x 5*^) respectively L2(M^) are given by 

{^gHx,R)U)(y,n) = U{g-' • (j,fi)) = U{R-'(y - x), R-'n), x,y e R^n e S^U e UiR^ x 5^), 
{iigHx,R)f)(y) = f{R-'(y-x)) , R e 50(3), x,j g R^/ g 



Intuitively, iig=(x,R) represents a rigid motion operator on images, whereas 2g=(^x,R) represents a rigid motion on 
HARDI images. 



Operator on HARDI-images must be left-invariant as the net operator on a HARDI-image should commute with 
rotations and translations. For detailed motivation see 1341 , where our motivation is similar as in our framework of 
invertible orientation scores I 



Tiieorem 1 Let JC be a bounded operator from L2(M^ x 5*^) into Loo(M^ x S'^) then there exists an integrable kernel 
k : {R^ x S'^) X (R^ x S'^) C such that ||/C|p = sup Jj^3y^s^ \k{y,n;y',n')\'^dyda{n') and we have 

{JCU){y,n)= I k{y,n',y',n')U{y',n')dy'dcj{n'), (19) 

for almost every iy^n) G x 5^ and all U G L2(M^ x S'^). Now JCk '= JC is left-invariant iffk is left-invariant, i.e. 

2goJCk=JCko2g<=>\/geSEi3)'^yyeR^^n,n'es^ ' ' (y ^n) ] Q - {y' ,n')) = k(y ,n ] y' ,n') . (20) 

Then to each positive left-invariant kernel /c : x 5^ x x 5^ ^ R+ with j^^ k{0^ez ; yjn)da{n)dy = Iwe 
can associate a unique probability density p :R^ y\ S'^ ^ R+ with the invariance property 

p{y,n) = p{Re^^ay^Re,,an), for all a G [0,27r), (21) 

by means ofp(y^ n) = /c(j, n ; Cz)- The convolution now reads 

JCkU{y,n) = {p^^s^s^U)(y,n) = J j p{Rl{y - y'),Rln)U{y' ,n')dcj{n')dy' , (22) 
where a denotes the surface measure on the sphere and where Rn' is any rotation such that n' = Rn'^z- 



14 



For details see 1341 . By the invariance property (21 ), the convolution (22) on x S'^ may be written as a (full) 
(3) -convolution. To this end we extend our positively valued functions U defined on the quotient >^ 6*^ = 
(R^ X SO{3))/{{0} X S0{2)) to the full Euclidean motion group by means of Eq. ^ which yields 

[/(X,i?e.,7^e„/3i^e.,a) = t/(x, n(/3, 7)) . (23) 

in Euler angles. Throughout this article we will use this natural extension to the full group. 

Definition 3 We will call : R^ x SO (3) R, given by Eq. pi, the HARDI- orientation score corresponding to 
HARDI image /7 : R^ x 5'^ ^ R. 

Remark 1 By the construction of a HARDI- orientation score, Eq. it satisfies the following invariance property 
U(x,RRe,^a) = U{x,R)forallx G R^i? G 5(9(3), a G [0,27r). 

An SE{3) convolution L18J of two functions p : SE{3) -^R,U : SE{3) ^ R is given by: 

{p^SEi3)U){9)= [ p{h-'g)U{h)d^sEi3){h)^ (24) 

JSE{3) 

where Haar-measure djj.sE{3){^^ R) = dxd/i5o(3)(^) with d/i5o(3)(^e;,,7^e^,/3^e^,cK) = sin /3dad/3d7. It is easily 
verified that that the following identity holds: 

(P*5£;(3) = 27r (p*R3>,52 U){x,Re^) . 

Later on in this article (in Subsection |8.1| and Subsection |8.2| ) we will relate scale spaces on HARDI data and first order 
Tikhonov regularization on HARDI data to Markov processes. But in order to provide a road map of how the group- 
convolutions will appear in the more technical remainder of this article we provide some preliminary explanations on 
probabilistic interpretation of R^ x 6*^ -convolutions. 

4 Left-invariant Vector Fields on SE{3) and their Dual Elements 

We will use the following basis for the tangent space Te{SE{3)) at the unity element e = (0, /) G ^£^(3): 

Ai d^jc.M = dy.A^, = dz.A^ d^.As O^.Aq = da , (25) 

where we stress that at the unity element (0, = /), we have (3 = and here the tangent vectors 9/5 and are not 
defined, which requires a description of the tangent vectors on the 5*0 (3) -part by means of the second chart. 

The tangent space at the unity element is a 6D Lie algebra equipped with Lie bracket 

[A,B] =lm^"^ {a{t)b{t){a{t))-\b{t))-^ - e) , (26) 

where t ^ a{t) resp. t ^ b{t) are any smooth curves in ^£^(3) with a(0) = 6(0) = e and a\0) = A and b\0) = B, 
for explanation on the formula ([26]) which holds for general matrix Lie groups, see ||29l App.G]. The Lie-brackets of 
the basis given in Eq. ([25]) are given by 



[A„A,] = ^4A, , (27) 

k=l 

where the non-zero structure constants for all three isomorphic Lie-algebras are given by 



^k^^k^i sgnperm{i - 3, j - 3,fe - 3} if ij,k > 4,i j k, 

I sgnpermjz, j — 3, k} if z, /c < 3, j > 4, z 7^ j 7^ k. 
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More explicitly, we have the following table of Lie-brackets: 





( 











Az 


-A2\ 


















Ai 


1,...6 











A2 


-Ai 





:1,...,6 ~ 





A3 


-A2 





Ae 


-As 









Ai 


-As 





A4 




V A2 


-Al 





As 


-A4 


/ 



(29) 



so for example cfg = 1, 0^4 = cfg = 0, c 



— — = — 1. The corresponding left-invariant vector fields {^^}^^^ 
are obtained by the push-forward of the left-multiplication Lgh = ghhy Ai\g(j) = {Lg)^Ai(j) = o L^) (for 

all smooth (j) : ^ R which are locally defined on some neighborhood of g) and they can be obtained by the 
derivative of the right-regular representation: 



A.\J={dn{A,)cj>){g) = ]im^ 



with ng(l){h) = (t){hg). 



(30) 



Expressed in the first coordinate chart, Eq. ([12]), this renders for the left-invariant derivatives at position 
g = (x,y,z,i?e^,^i?ey,/3^e,,a) ^ SE{?>) (scc also |18, Scction 9.10]) 



Aa 

Aq 



(cos a cos P cos 7 — sin a sin 7) dx + (sin a cos 7 + cos a cos /3 sin 7) 



(— sin a cos /3 cos 7 — cos a sin 7) dx + 

sin /3 cos dx + 

cos a cot P da + 

— sin a cot /3 + 



(cos a cos 7 — sin a cos /3 sin 7) dy 
sin /3 sin 7 
sin a 
cos a d B 



— cos a sin /3 , 
+ sin a sin /3 , 
+ cos dz, 

cos Oi f) 



+ 



sin /3 
sin q; o 
sin ^ ^7 5 



for /3 7^ and /3 7^ tt. The explicit formulae of the left-invariant vector fields in the second chart, Eq. ( 15 ), are : 

Ai — cos a cos P dx + (cos 7 sin a + cos a sin /3 sin 7) dy Aa — — cos o; tan f3 da -\- sin a + ^7 , 

+ (sinasin7 — cos a cos 7 sin ^5) dz, 
A2 — — sin a cos (3 dx + (cos a cos 7 — sin a sin f3 sin 7) 

+ (sin a sin ^5 cos 7 + cos a sin 7) dz , 
^3 = sin P dx — cos /3 sin j dy -\- cos /3 cos 7 , Ae = da 



A5 — sin a tan /3 da -\- cos a ^5 



P cos^ 



sm q; o 
5 



(31) 



(32) 



for ;5 7^ I and f3 Note that d7^ is a Lie-algebra isomorphism, i.e. 



[Ai.Aj] = ^c^jAk ^ [d7^(A,),d7^(A,)] = ^4•d7^(Afe) ^ [A^Aj] = AiAj - AjAi = ^4^^ 



k=l 



k=l 



k=l 



These vector fields form a local moving coordinate frame of reference on SE{3). The corresponding dual frame 
{d^^ , . . . , d^^} G {T{SE{3))y is defined by duality. A brief computation yields : 

/ d^M / \ 

d^' I / ^ ^ XT . X I d2/ 

dz 
da 



dA^ 
dA^ 
d^^ 
V ^A' J 














(Re^^^R^^jRe^^a) 


^ 








V d7 / 

where the 3 x 3-zero matrix is denoted by and where the 3 x 3-matrices Mp^^^ - are given by 



( dx\ 
dy 
dz 
da 

dp 

V d7 y 



(33) 




-cosasin/3 \ / -cosatan/3 sina 

sin a sin /3 , M^^^ = sin a tan /3 cos a 
cos/3 / V 1 




Finally, we note that by linearity the z-th dual vector filters out the i-th component of a vector field Yl^j=i '^^Aj 



{dA\^v^Aj 



for alH = 1, . . . , 6. 
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Remark 2 In our numerical schemes, we do not use the formulas {31) and {32) for the left-invariant vector fields 



as we want to avoid sampling around the inevitable singularities that arise with the the coordinate charts, given by 



Eq. {12) and ([75]), of S^. Instead in our numerics we use the approach that will be described in Section 12 However, 
we shall use formulas p7] ) and p2]) frequently in our analysis and derivation of Green's functions of left-invariant 
diffusions and left-invariant Hamilton-Jacobi equations on y\ S'^. These left-invariant diffusions and left-invariant 
erosions are similar to diffusions and erosions on R^, we ''only'' have to replace the fixed left-invariant vector fields 
, . . . , } by the left-invariant vector fields { |^ , . . . , ^5 1^} which serve as a moving frame of reference (along 
fiber fragments) in R^ x S*^ = SE{?>) / {{0} x S0{2)). For an a priori geometric intuition behind our left-invariant 
erosions and diffusions expressed in the left-invariant vector fields see Figure^ 



5 Morphological Convolutions on x 

Dilations on the joint space of positions and orientations R^ x S'^ are obtained by replacing the (+, - j-algebra by the 



(max, +)-algebra in the R^ x 6'^ -convolutions (22) 



(fc©M3x5^ /7)(y,n) = sup [fc(i?J(y -/), i^Jn) + n^] (34) 

(y',n')GM3xS'2 

where k denotes a morphological kernel. If this morphological kernel is induced by a semigroup (or evolution) then 
we write kt for the kernel at time t. Our aim is to derive suitable morphology kernels such that 

h ©R3x52 ks ks+t , for all 5, t > 0, 

where t ^ /c^(y, n) describes the growth of balls in R^ x 6*^, i.e. t ^ kt{y^n) is the unique viscosity solution, see 
1201 and Appendix|B] of 

with the inverse of the metric tensor restricted to the sub-Riemannian manifold {SE{3)^ d^^, d^^) 

G = -^(^1 (g) ^1 + A (g) A2) + -^(^4 g) ^4 + A g) A), (36) 

given by 

G"^ = D'^^{Ai (g) ^1 + A (g) A2) + D^^{A4 g) ^4 + A A)- 
Both G and its inverse are well-defined on the cosets R^ x 5^ = (R^ x 6'O(3))/({0} x S0{2)), cf. Appendix|El 



Furthermore, in Eq. (35 ) we have used the morphological delta distribution S^{y^ n) = +00 if (y, n) ^ (0, e^) and 



HmA:^ ©r3>,52 U = {-S^) ©r3>,52 U = U 



else, where we note that 

uniformly on R^ x S*^ . Furthermore, we have used the left-invariant gradient which is the co-vector field 

6 5 

dC/(y,n) = ^(A(C/))(y,n) d^'|(^_„^ = ^(A(f/))(y, n) d^^l^^^^ , (y,n) eR" ^ S', (37) 

i=l ' i=l 

which we occasionally represent by a row vector given by 

VC/(y, n) = (AiU iy,n),..., A^U{y, n) , 0) . (38) 
The dilation equation (|35]> now becomes 



f^(y,n,t) = i {D'HiAiW{y,n,t)f + {A2W{y,u,t))^) + D^HiAiW{y,u,t)f + {A5W{y,n,t)f)) 
W{y,n,0) =-5^{y,n). 
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Now we can consider either a positive definite metric (the case of dilations), or we can consider a negative definite 
metric (the case of erosions). In the non-adaptive case this means; either we consider the gu = > or we choose 
them gu = < 0. Note that the erosion kernel follows from the dilation kernel by negation kt ^ —kt. Dilation 

kernels are negative and erosion kernels are positive and therefore we write 

kt := -h > and k^ := kt < . 

for the erosion kernels. Erosions on x S'^ are given by: 



(fc+eR3,s^[/)(y,n)= inf [[/(y',n') + (y - y'),i??n)] . 

We distinguish between three types of erosions/dilations on x S'^ 



(39) 



1. Angular erosion/dilation (i.e. erosion on glyphs): In case gn = gss = and ^44 < the erosion kernels are 
given by 

ify/0 



i^t iy^nj - I ^+,^11=0, 

such that the angular erosions are given by 



1^441 



(0, n) else 



= inf 



{0,R^,in')) + U{y',n') 



whereas the angular dilations ^44 > are given by 



sup 



^^1^1=^33=0,^44 (Q^^T^) + U{y,n') 



(40) 



(41) 



2. Spatial erosion/dilation (i.e. the same spatial erosion for all orientations): In case gugss ^ 0, gu < 0, gss < 
and ^44 = the erosion kernels are given by 



^ + ,5^11, 5^33, 5^44=0 



= I ^,l^ll| = l^33h 



^44=0 



if n 7^ 



(y,e;,) else 



such that the spatial erosions are given by 



/7,+ ,|fi'll|, 1^331, ^44=0 ^ TT\f^T n\ — /iL + 'l^lll'l^33|,^44=0/ „ \ ^ TTf nWf^r^ 

whereas the spatial dilations gu^gss > 0, gugss > are given by 



,911,933,944=0 



■.^s- U){y,n) = (fcr'^"'^^^'^^^=°(-,e,) 0^3 f/(-,n))(y). 



3. Simultaneous spatial and angular erosion/dilations (i.e. erosions and dilations along fibers). The case ^44 7^ 
and ^11 7^ or ^33 ^ 0. 



Similar to our previous work on R^ y\ /S^ -diffusion f26l the third case is the most interesting one, simply because one 
would like to erode orthogonal to the fibers such that both the angular distribution and the spatial distribution of the a 
priori probability density U : R^ x 5^ ^ R+ are sharpened. See Figure [s] 
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Figure 8: Top row: output angular erosion on a single glyph with crossings for several parameter settings. Middle row: 
input and output simultaneous spatial and angular erosion (i.e. x 6*^ -erosion), where angular erosion dominates 
t = 3, D^^ = 0.6, D^^ = 0.3. Bottom row: input and output >^ iS^-erosion, where angular spatial erosion dominates 
t = 2,D^^ = 0.1, = 0.3. 
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6 Left-Invariant Diffusions on SE{S) = R^ SO{3) and x 



In order to apply our general theory on diffusions on Lie groups, 1271 . to suitable (convection-)diffusions on HARDI 
images, we first extend all functions U : xi S'^ ^ R+ to functions ^/ : x 5*0 (3) R+ in the natural way, by 
means of ([3]). 

Then we follow our general construction of scale space representations W of functions U (could be an image, or a 
score/wavelet transform of an image) defined on Lie groups, |27|, where we consider the special case G = SE{3)'. 



r dtW{g, t) = Q^'^(^i, ^2, . . . , A) W{g, t) , 
I \^W{g,t) = U{g). 

which is generated by a quadratic form on the left-invariant vector fields: 

6 6 

i=i j=i 



(42) 



(43) 



Now the Hormander requirement, BSl . on the symmetric D = [D'^^] G R^^^, D > and a, which guarantees smooth 
non-singular scale spaces, for SE{3) tells us that D need not be strictly positive definite. The Hormander requirement 
is that all included generators together with their commutators should span the full tangent space. To this end for 
diagonal D one should consider the set 

<S = {iG{l,...,6} |I)^VOVa,^0}, 

now if for example 1 is not in here then 3 and 5 must be in S, or if 4 is not in S then 5 and 6 should be in 5. 
If the Hormander condition is satisfied the solutions of the linear diffusions (i.e. D, a are constant) are given by 
(3) -convolution with a smooth probability kernel p^'^ : SE{3) R+ such that 

W{g,t) = (pf ^sEis) U){g) = I pf^\h-^g)U{h)d^sEi3){h), 

SEiS) 

UO ^^^'^ *'^^^^^ U = U , withpf'^ > and J^^^^^p^'^{g)d^sEi3){g) = 1- 

where the limit is taken in L2(5'£^(3))-sense. On HARDI images whose domain equals the homogeneous space 
R^ y\ S'^ one has the following scale space representations: 



dtW{y, n, t) = Q^(U)MU) (^,, • • • , A) W{y, n, t) , 
I^(y,n,0) = /7(y,n). 



(44) 



withQ^(^)'«(^)(A, A,..., A) = ELi (^iA + E^=iA^'^(^) a), where from now on we assume that D(/7) 
and a([/) satisfy 

a{U){gh) = Z^{a){g) and D(/7)(^/i) = Z^D{U){g)Z^ . (45) 
for all g G SE{3) and all h = (0, i?e,,a) ^ ({0} x S0{2)) and where 
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— sin a 
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1 ) 



(46) 



Recall the grey tangent planes in Figure[6] where we must require isotropy due to our embedding of R^ x S'^ in 5'^(3), 
cf. m ch:4]. 
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In the linear case the solutions of (44) are given by the following kernel operators on x: 5^: 



W^(y,n,t) = (pf''' *R3>,s. U){y,n) 

= l7l pf'''((i?e„yi?e„/3')^(y - y'), (i?e.,yi?e„/3')^n)) C7(y', n(/?', 7' )) dy'dcT(n(/?',7')), 

R3 



where the surface measure on the sphere is given by dcr(n(/3', 7')) = sin f3^ d^d^^ = da{h{$^ t)) = I cos $\ d;^d7. 
Now in particular in the Hnear case, since (R^, /) and (0, SO (3)) are subgroups of SE{3), we obtain the Laplace- 
Beltrami operators on these subgroups by means of: 

As2 = Q^=diag{0,0,0,l,l,l},a=0 ^ (^^)2 ^ (^^)2 ^ (^^)2 ^ (5^)2 ^ cot(/3)a^ + siR-^ (9^)^ , 
Am3 = Q^=d-g{l,l,l,0,0,0},a=0 ^ (^^)2 ^ (^^)2 ^ (^3)2 ^ (5^)2 ^ (5^)2 ^ (5^)2 ^ 

One wants to include line-models which exploit a natural coupling between position and orientation. Such a coupling 
is naturally included in a smooth way as long as the Hormander's condition is satisfied. Therefore we will consider 
more elaborate simple left-invariant convection, diffusions on SE{3) with natural coupling between position and 
orientation. To explain what we mean with natural coupling we shall need the next definitions. 



Definition 4 A curve 7 : R+ ^ x 6*^ given by s ^ 7(5) = {y{s)^n{s)) is called horizontal if n{s) = 
\\y{s)\\~^y{s). A tangent vector to a horizontal curve is called a horizontal tangent vector A vector field A on 
R^ XI S'^ is horizontal if for all (y^n) G R^ x S'^ the tangent vector A{y^n) horizontal. The horizontal part Hg of 
each tangent space is the vector- sub space ofTg{SE{3)) consisting of horizontal vector fields. Horizontal diffusion is 
diffusion which only takes place along horizontal curves. 



It is not difficult to see that the horizontal part Hg of each tangent space Tg{SE{3)) is spanned by {^3, ^4, ^5}. 
So all horizontal left-invariant convection diffusions are given by Eq. ( [44| ) where in the linear case one must set 
ai = a2 = aQ = 0, Dj2 = D2j = Dij = Dji = Djq = Dqj = for all j = 1, 2, . . . , 6. Now on a commutative 
group like R^ with commutative Lie-algebra {S^^^, . . . ^Oxq} omitting 3-directions (say dx^^, dx^, Oxq) from each 
tangent space in the diffusion would yield no smoothing along the global xi, X2, xe-axes. In SE{3) it is different 
since the commutators take care of indirect smoothing in the omitted directions {.Ai, .42, ^e}, since 

span {^3,^4,^5, [^3, A] = ^2, [A, A] = A, [As, A3] = Ai} = T{SE{3)) 

Consider for example the iSE^ (3) -analogues of the Forward-Kolmogorov (or Fokker-Planck) equations of the direction 
process for contour-completion and the stochastic process for contour enhancement which we considered in our pre- 
vious work, (301 , on SE{2). Here we first provide the resulting PDEs and explain the underlying stochastic processes 
later in subsection|8.1[ The Fokker-Planck equation for (horizontal) contour completion on SE{3) is given by 



dtWiy,n,t) 
limH^(y,n,t) 



(-^3 + D{{Ai)^ + (A)')) W{y, n, t) = {-A3 + D A<j2) W{y, n,t) ,D 
U{y,n). 



> 0. 



where we note that (^6)^(^(y, n(/3, 7), s)) = 0. This equation arises from Eq. (44) by setting D^^ = D^^ 
as = 1 and all other parameters to zero. The Fokker-Planck equation for (horizontal) contour enhancement is 

r dtW{y,n,t) = {D^^iAs)' ^ D^\{A,)'^{A,)'))W{y,n,t) = {D^^iAs)' ^ D As^) W{y,n,t) , 
j \imW{y,n,t) =U{y,n). 



(48) 
i:>and 



(49) 



The solutions of the left-invariant diffusions on R'^ x S given by ( 44 ) (with in particular ( 48 ) and ( 49 )) are again 



given by convolution product (147]) with a probability kernel on ! 
kernels, see Figure fTOl 



y\ S . For a visualization of these probability 
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7 Left-invariant Hamilton- Jacobi Equations on x 5^ 



The unique viscosity solutions of 

\ l^(y,n,0) = C/(y,n), 



where r] > |, are given by dilation, Eq. ( |34[ ), with the morphological Green's function k 

M^(y,n,i) = (fcf"'^"'''-©f/)(y,n), 
whereas the unique viscosity solutions of 



t 



(y,n,t) = -2^G^yi„) (dl^(-, •,t)|y^„ , dM/(-, •,t)|y,„))'' 
Vr(y,n,0) = C/(y,n) 



at 



(50) 



(51) 



(52) 



are given by erosion, Eq. (39 ), with the morphological Green's function '^'^ : x — > 

H^(y,n,t) = (fcf"^"'''+eC/)(y,n), 

This is formally shown in Appendix [BJ Theorem]?] The exact morphological Green's functions are given by (where 
we recall ( flOl l) 



(53) 



-">-"--(y,„) = fcf">-"-+(y,„) := 



mf / £„(7(p),7(p)) f ^ I dp. 

7(0) = (0,1 = Re^),7(t) = (y, i^n), Q 



(54) 



with Lagrangian 



'Cr,(7(p),7(p)) := 



2r]- 1 



where we applied short notation 7*(p) = (d^* |^^^^ , j{p)) and with R"^ y\ 5 -"erosion arclength" given by 



P{r) =/VG^(f)(7(r),7(r))df = / / E ^Kd^^l^(f) , 7(r))P dr. (55) 

Y ^^{1,2,4,5} 

For further explanation an details see Appendix |B] in particular Lemma [2] As motivated in Appendix [B] we use the 
following asymptotical analytical formula for the Green's function 



(y,ii(/3,7)) = ± 



2r] 



|P(y,a = Q,^,7)P 



(56) 



for sufficiently small time t > 0, where C > 1 is a constant that we usually set equal to 1 and where the constants 



~ ^, i = 1, . . . , 6, are components of the logarithm 



that we shall derive in Section [971] Eq. (75 ). 



Apparently, by Eq. (54) the morphological kernel describes the growth of "erosion balls" in x 5^. In Appendix 
|B]we show that these erosion balls are locally equivalent to a weighted modulus on the Lie-algebra of SE{3), which 
explains our asymptotical formula (56). This gives us a simple analytic approximation formula for balls in x] 6*^, 



where we do not need/use the minimizing curves (i.e. geodesies) in (54 ). 
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adaptive 
erosion Su. 
dilation 



" • • • 



c=-2 c = -l c=0 c=l c=2 




Figure 9: From left to right, input glyph at say y = 0. Sample points where Alb^(0, n) < resp. > are respectively 
indicated in red and blue. Linearly eroded glyph {D^^ = 0.4, t = 0.5, r] = 1, At = 0.01), adaptive eroded glyph 
(D^^ = 0.4, t = 0.5, r] = h (Pix) = D^^x^). 



Remark 3 In Appendix\G\we provide the system ofPfaffian equations for geodesies on the sub-Riemannian manifold 
{SE{3)^ dA^^ dA'^j dA^) as a first step to generalize our results on SE{2) = x 6'^ in [3L App. A], where we 
generalize our results on SE{2) = M? xi in 07] App. A] to {SE{3). In order to compute the geodesies on the 
sub-Riemannian manifold {SE{3),dA^^ dA^) (used in the erosions Eq. 52) a similar approach can be followed. 



7.1 Data adaptive angular erosion and dilation 



In the erosion evolution (52 ) one can include adaptivity by making I)^^ depend on the local Laplace-Beltrami-operator 

D^\U){y,n) = 4>{ALBU{y,n)-c). 

with (j) a non-decreasing, odd function and inf Ai,^/7(y, n) < c < sup Ai^^[/(y, n). The intuitive idea 

here is to dilate on points on a glyph x + /i{[/(x,x)n \ n e S'^} where the Laplace-Beltrami operator is negative 
(usually around spherical maxima) and to erode on at locations on the glyph where the Laplace-Beltami operator is 
positive. Parameter c > tunes the boundary on the glyph where the switch between erosion and dilation takes place. 
See Figure [9] where we have set = I)^^ while varying c. 



8 Probability Theory on x 

8.1 Brownian Motions on SE{3) = >^ S0{3) and on >^ 

Next we formulate a left-invariant discrete Brownian motion on SE{3) (expressed in the moving frame of reference). 
The left-invariant vector fields {^1, ... ,^6} form a moving frame of reference to the group. Here we note that 
there are two ways of considering vector fields. Either one considers them as differential operators on smooth locally 
defined functions, or one considers them as tangent vectors to equivalent classes of curves. These two viewpoints are 
equivalent, for formal proof see I?] Prop. 2.4]. Throughout this article we mainly use the first way of considering vector 
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fields, but in this section we prefer to use the second way. We will write {ei (^) , . . . , ee (^) } for the left-invariant vector 
fields (as tangent vectors to equivalence classes of curves) rather than the differential operators {^i|^ , • • • , ^el^}- 
We obtain the tangent vector from Ai by replacing 



^ (1,0,0,0,0,0), 
dy ^ (0,1,0,0,0,0), 
d, ^ (0,0,1,0,0,0), 



d(3 ^ (0, 0, 0, a cos P cos 7, a cos P sin 7, —a sin /3), 

^ (0,0,0, a cos 7, a sin 7,0), 
da ^ (0, 0, 0, cos7sin/3, sin 7 sin /3, cos /3), 



(57) 



where we identified 5*0 (3) with a ball with radius 27r whose outer-sphere is identified with the origin, using Euler 
angles i?e^,7^e^,/3^e^,Q! ^ cen(/3, 7) G 5o,27r. Ncxt wc formulate l eft-in variant discrete random walks on SE{3) 
expressed in the moving frame of reference {e^}^^^ given by ( 31 ) and (112): 



(Yn+l,Nn+l) = (Yn,Nn) + As 
(Yo,no)-/7^, 



5 5 

2=1 j=l 



foralln = 0, ...,A^-1, 



with random variable (Yq, Hq) is distributed by U^, where are the discretely sampled HARDI data (equidistant 
sampling in position and second order tessalation of the sphere) and where the random variables (Y^,Nn) are re- 
cursively determined using the independently normally distributed random variables {^i,n+i}r=r/'... 5^~^' ^i,n+i ~ 
A/'(0, 1) and where the stepsize equals As = and where a := Ei=i ^i^i denotes an apriori spatial velocity vector 



having constant coefficients with respect to the moving frame of reference {e^lf^^ (just like in (43)). Now if we 
apply recursion and let ^ oo we get the following continuous Brownian motion processes on SE{3)'. 



Y{t) = Y{0) + / ( E e 

\i=l 



N{t) = N{0) + / ( X: a, e, 

Vi=4 



i|(y(r),Ar(r)) " 
(y(r),iV(r)) 



5 



r 2 



J=4 



-Jl(y(r),iV(r)) 



dr. 



(58) 



with Ei - A/'(0, 1) and (X(0), N{0)) - U and where a = V2D e 



3x6 



cr > 0. Note that d^/r ^ 



"2 dr. 



Now if we set U = ^o,e^ (i-e. at time zero ) then suitable averaging of infinitely many random walks of this process 
yields the transition probability (y, n) p^'^(y, n) which is the Green's function of the left-invariant evolution 
equation s ([44| ) on x S'^. In general the PDE's (44 ) are the Forward Kolmogorov equation of the general stochastic 
process ( |58l). This follows by Ito-calculus and in particular Ito's formula for formulas on a stochastic process, for 
details seelll app.A] where one should consistently replace the left-invariant vector fields of by the left-invariant 
vector fields on R^ x 5^. 



In particular we have now formulated the direction process for contour completion in R^ y\ S'^ (i.e. non-zero parameters 
in (58) are I)^^ = D^^ > 0,a3 > with Fokker-Planck equation (48)) and the (horizontal) Brownian motion for 
contour- enhancement in R^ x 5*^ (i.e. non-zero parameters in (58 ) are D^^ > 0, D^^ = D^^ > with Fokker-Planck 



equation (49)). 



See Figure [T0| for a visualization of typical Green's functions of contour completion and contour enhancement in 
R^ X 6'^-\d = 2,3. 



8.2 Time Integrated Brownian Motions 

In the previous subsection we have formulated the Brownian-motions ( [58] ) underlying all linear left-invariant convection- 
diffusion equations on HARDI data, with in particular the direction process for contour completion and (horizontal) 
Brownian motion for contour-enhancement. However, we only considered the time dependent stochastic processes 
and as mentioned before in Markov-processes traveling time is memoryless and thereby negatively exponentially 
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R2 X, s- 



Contour completion 





R3 X 



^3—1 



I<C?1 



■my 




i^33 = l 

/?44 ^i^55>0 



Contour enhancement 
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Figure 10: The Green's function of Forward Kolmogorov equation of (the horizontal) direction process and the Green's 
function of the Forward Kolmogorov equation of (horizontal) Brownian motion in y\ S^~^, for d = 3 as considered 
in this article and for d = 2 as considered in our previous work 1331 (completion) and ||30l[3T| (enhancement). 
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distributed T - NE{X), i.e. P(T = t) = Xe'^^ with expectation E(T) = A-\ for some A > 0. By means 
of Laplace-transform with respect to time we relate the time-dependent Fokker-Planck equations to their resolvent 
equations, as at least formally we have 

W{y,n,t) = (e^(^'"(^))[/)(y,n) andP^(y,n,t) = A / e-^^(e^(^'"^^))/7)(y, n) = A(A/ - Q'^'^(V))-i[/(y, n), 

Jo 

for A > and all y G M^, n G 5*^, where the negative definite generator Q^'^ is given by (43) and again with 
\/U = {AiU, . . . , AgU)^. The resolvent operator A(A/ - Q^=diag(i)^^),a=0(^y ^ ^-i occurs in a first order Tikhonov 
regularization as we show in the next theorem. 

Theorem 2 Let U G L2(M^ x S'^) and A, D^^ > 0, D^^ = D^^ > 0. Then the unique solution of the variational 
problem 

argmin / ^(P(j,n) - U{y,n)f + y D^^\AkP{y,n)\^dyd(j{n) (59) 

is given by P^{y^n) = {R^ *r3x52 U){y^n), where the Green's function : y\ S'^ ^ is the Laplace- 

oo 

transform of the heat-kernel with respect to time: R^ (j, n) = X J (j, n)e~^^ dt with D = diag{D^^^ . . . , D^^}. 



P^{y^n) equals the probability of finding a random walker in y\ S'^ regardless its traveling time at position y G 
with orientation n G S'^ starting from initial distribution U at time t = 0. 

For a proof see |[26ll . Basically, P^ (y, n) = {R^ *r3 /7) (y, n) represents the probability density of finding a random 
walker at position y with orientation n given that it started from the initial distribution U regardless its traveling time, 
under the assumption that traveling time is memory less and thereby negatively exponentially distributed T ~ NE{X). 
There is however, a practical drawback due to the latter assumption: Both the time-integrated resolvent kernel of the 
direction process and the time-integrated resolvent kernel of the enhancement process suffer from a serious singularity 
at the unity element. In fact by some asymptotics one has 

^„^^a=(0,0,l),D=diag{0,0,0,D«,D-}(Q^Q^^^^^ ^ ^ for < Z « 1, 

where \g\D^^,D^^,D^^ is the weighted modulus on SE{3) = x 50(3). For details see Appendix [d| These kernels 
can not be sampled using an ordinary mid-point rule. But even if the kernels are analytically integrated spatially and 
then numerically differentiated the kernels are too much concentrated around the singularity for visually appealing 
results. 



8.2.1 A /c-step Approach: Temporal Gamma Distributions and the Iteration of Resolvents 



The sum T of k independent negatively exponentially distributed random variables Ti ^ NE{X) (all with expectation 
E{Ti) = A"^) is Gamma distributed: 

k \ k -i-k X 

T = ^ T, with pdf : p{T = t)= p(Ti = t) p(T, =t) = T(t ; A) := Jj^^f'^' , ^ ^ N, 

where we recall that temporal convolutions are given by (/ g){t) = Jq f(t — r)g{r) dr and note that application 
of the laplace transform C, given by Cf{X) = f{t)e~^^dt yields C{f ^]^+ g) = C{f) • C{g) . Now for the sake 
of illustration we set /c = 2 for the moment and we compute the probability density of finding a random walker with 
traveling time T = Ti + T2 at position y with orientation n given that it at t = started at (0) with orientation 
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e^. Basicly this means that the path of the random walker has two stages, first one with time Ti ~ NE{X) and 
subsequently one with traveling time T2 ~ NE{X) and straightforward computations yield 

00 

<t=2(y. n) = / p{Gt = (y, n)\T = t and Go = e) p{T = t) dt 



I p{Gt = {y,n)\T = Ti+T2 = t and Go = e) p(Ti +T2 = t) dt 


00 t 

J J p{Gt,+T2 = {y,n)\Ti=t-s and T2 = s and Go = e) p(ri = t - s) p(r2 = s) dsdt 







A2 £ t ^ /(i^i 



5.2 *R3x,s2 (5e)(y,n)ds (A) 



A2 £ ( t ^ J{Kt-s *K3>,52 ii:,)(y, n)ds ) (A) 



X'C{t^ Kt{-)) (A) *R3>,s^ £(i ^ (A)(y,n) = (i?^t=i *R^>.s^ Kl=i)iy^«) ■ 



pD,a 



(60) 



By induction this can easily be generalized to the general case where we have 

<'L2 = ^a'" *m^^5^ ^a'" With i?^'" = Rll^, mdp{T = t) = {p{T^ = •) p{n = -W) ■ 

As an alternative to our probabilistic derivation one has the following derivation (which holds in distributional sense): 

00 

R^^'l = J{e'Q''''(^)5e)r{t;k,X)dt 

= ^(-Q°'''(V ) + A/)-'=5e = (A {-Q'^'^iiV ) + M)-')''Se 

_ pD,a k — 1 pD,a 

where we note that the Laplace transformation of a Gamma distribution equals 

£r(-,fc,A)(5) = + 



(A,fc) = (l,4) 




Figure 11: Glyph visualization, recall Def. [ij of the kernels i?^'^^2 • >^ ^S*^ with diffusion matrix 

D = diagjO, 0, 0, 1)^^, I)^^, 0} and convection vector a = (0, 0, 1, 0, 0, 0), for several parameter-settings (A, k) for 
1)44 _ Q QQ5 Kernels are sampled and computed on a spatial 8 x 8 x 8-grid and on an a 162-point tessellation of the 
icosahedron using a left-invariant finite difference scheme, cf. | 21 1. For the sake of comparison we fixed the expected 
value of the Gamma-distribution to | = 4. The glyph- visualization parameter /i (which determines the global scaling 
of the glyphs) has been set manually in all cases. 
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(A,fc) = (i 1) 



(A,fc) = (1,2) 



input 
(A,fc) = (1,0) 





(A,fc) = (§,3) 





(A,fc) = (2,4) 




Figure 12: Glyph visualization, recall Def. [T] of left-invariant finite difference scheme, onal2xl2x 12-grid 
and 162-point tessellation of the icosahedron, applied to the input data set, according to our time integrated contour- 
completion process, ref. Eq. (58 ), with diffusion matrix D = diag{0, 0, 0, D^^, D^^}, D^^ = 0.005 and convection 



vector a = (0, 0, 1, 0, 0), iterated A;-times with E{Ti) = A"\ i = 1, 



Figure 12 shows some experiments of contour completion on an artificial data set containing fibers with a gap we 
would like to complete, for various parameter settings of (A,/c), where Ti ~ NE(A) and k denotes the iteration 
index of the time integrated contour-completion process. In principle this scheme boils down to x 5^ -convolution 
with the Green's functions depicted in Fig. [13] (this is only approximately the case due to discretization). For a fair 
comparison we kept the expected value constant, that is £^(T) = £^(Ti + . . . T^) = | = 4 in Figure [ll] and E{T) = 



E(Ti^...n) 



2 in Figure 
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which roughly coincides with half the size of the gap. Note that the results 
hardly change after two iterations, as the graphs of the scaled Gamma distributions r(- ; |, k) /V{{2(k — l)/k)] |, /c) 
are similar for /c = 2, 3, 4, 5. 



8.3 Cost Processes on SE{?>) 

In this subsection we present a short overview of cost processes on SE{?>). The mapping / ^ lim^^o loge ^ defines a 
morphism of the (+, •) -algebra to the (min, +) -algebra on R+ (the co-domain of our HARDI orientation scores, recall 
Def. [s]) and it is indeed readily verified that lim^^o log^ • = a + 6 and lim^^o loge(^^ + — min(a, h). Using 
this map, various notions of probability calculus can be mapped to their counterparts in optimization problems. Next 
we mention the following definitions as given by Akian, Quadrat and Viot in adapted to our case. 



Definitions The decision space is the triplet {SE{3),0{SE{3)),C) with 0{SE{3)) denoting the set of all open 



subsets of the topological space SE{3). The function C : 0{SE{3)) 



U {oo} is such that 



1. C{SE{3)) = 
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2. C{4>) = +00 

3. C(U„ A„) = inf„ C{A„)forany An e 0{SE{3)) 

C is called a cost measure on SE{3). The map c : SE{3) given by g ^ c{g) such that C{A) = inf^^^ c{g) 

for all A C SE{3), is called the cost density of the cost measure C. 



Definition 6 Analogous to the random variables of probability theory, a decision variable G (on {SE{3), 0{SE{3))^ C) ) 
is a mapping from SE{3) to R. It induces a cost measure Co on R given by Cg{I) = C{G~^{I)) for all I G 0(R). 
The associated cost density is denoted by cq- 



One can formulate related concepts such as independent decision variables, conditional cost, mean of a decision vari- 
able, characteristic function of a decision variable etc. in the same way as in probability theory keeping in mind the 
morphism of the (+, •) -algebra to the (min, +) -algebra. The Laplace or Fourier transform in (+, •) -algebra corre- 
sponds to the Frenchel transform in the (min, +) -algebra. Now we present the decision counterpart of the Markov 
processes, namely Bellman processes. 



Definition 7 A continuous time Bellman process Gt on {SE{3)^ 0{SE{3)), C) is a function from SE(3) to C(R+) 
(set of continuous functions on non negative reals) with cost density 

POO 

ccili-)) = co(7(0)) + / c(t, 7W, f W) dt (61) 

where c is called the transition cost which is a map from R x SE{3) x T{SE{3)) to R such that 

inf c{t,g, A\ ) = for all {t,g) G R+ x 5'^(3), 

and where cq is some cost density on SE{3). 



We set 

c(t,7(t),7'(t))= J2 \{dA%^,yl'm\ 

iG{l,2,4,5} 

Then the marginal cost for a Bellman process Gt on {SE{3)j 0{SE{3))j C) to be in a state g at time t given an initial 
cost Co, rh{g^ t) := C{Gt = g) satisfies the following relation known as the Bellman equation 

dtfa + {Tc{SE{?,)){L)){dfa) = 0, 

rh{g,0)=co{g) ^^^^ 

where J^jC(se(3)) denotes the Fenchel transform, see Appendix [b| Definition [9j on the Lie- Algebra C{SE{3)) of 
left-invariant vector fields and where 

L{A)= E ^\{dA\A)\\ 

iG{l,2,4,5} 

dm = E Aim dA^ 

iG{l,2,4,5} 

for all A left-invariant vector fields within contact-manifold {SE{3)^ dA^^ dA^) and with left-invariant gradient of rh 
within {SE{3)^ dA^^ A^). This Bellman equation for the cost process coincides with the Hamilton- Jacobi equation on 
SE{3) ( |136| ) whose viscosity solution is given by morphological convolution with the corresponding morphological 
Green's function as proven in Appendix [B] 
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9 Differential Geometry: The underlying Cartan-Connection on SE{3) and 
the Auto-Parallels in SE{3) 



Now that we have constructed all left-invariant scale space representations on HARDI images, generated by means of 



a quadratic form (43 ) on the left-invariant vector fields on SE{3). The question rises what is the underlying differential 



geometry for these evolutions ? 

For example, as the left-invariant vector fields clearly vary per position in the group yielding a moving frame of 
reference attached to luminosity particles (random walkers in x S'^ embedded in SE{3)) with both a position and 
an orientation, the question rises along which trajectories in xi S'^ do these particles move ? Furthermore, as the 
left-invariant vector fields are obtained by the push-forward of the left-multiplication on the group, 

Ag = {Lg)^Ae, i.e. Ag4> = Ae{4> o Lg), where Lgh = gh, g,h e SE{3), ^ : SE(3) R smooth , 

the question rises whether this defines a connection between all tangent spaces, such that these trajectories are auto- 
parallel with respect to this connection ? Finally, we need a connection to rigid body mechanics described in a moving 
frame of reference, to get some physical intuition in the choice of the fundamental constant^ {aijf^^ and {^*^ }f^j=i 



within our generators (43 ). 



In order to get some first physical intuition on analysis and differential geometry along the moving frame {^i , . . . , 
and its dual frame {d^^, . . . , d^^}, we will make some preliminary remarks on the well-known theory of rigid body 
movements described in moving coordinate systems. Imagine a curve in R^ described in the moving frame of reference 
(embedded in the spatial part of the group SE{3)), describing a rigid body movement with constant spatial velocity 
c^^^ and constant angular velocity c^'^^ and parameterized by arc-length 5 > 0. Suppose the curve is given by 

3 

y{s) = ^a\s) A|y(,) where G C\%L],R), 

i=l 

such that ^ = Xli^i ^i\y{s) s > 0. Now if we differentiate twice with respect to the arc-length 

parameter and keep in mind that ^ ^ily(s) ~ ^ ^*ly(s)' 

y{s) = + 2c^^^ X c^^^ + c^^^ X (c^^^ x y{s)) . 

In words: The absolute acceleration equals the relative acceleration (which is zero, since c*^^^ is constant) plus the 
Coriolis acceleration 2c^ ^ x and the centrifugal acceleration x (c^ x y{s)). Now in case of uniform circular 
motion the speed is constant but the velocity is always tangent to the orbit of acceleration and the acceleration has 
constant magnitude and always points to the center of rotation. In this case, the total sum of Coriolis acceleration and 
centrifugal acceleration add up to the well-known centripetal acceleration, 

y{s) = 2c^'^ X (-c^') X Rr{s))^c^ x (c^'^ x Rr{s)) = -\\c^^^fRr{s) = -^^r(s), 

R 

where R is the radius of the circular orbit y{s) = m + i?r(s), ||r(5)|| = 1). The centripetal acceleration equals half 
the Coriolis acceleration, i.e. y{s) = c*^^^ x c^^\ 



In our previous work |[30l part II] on contour-enhancement and completion via left-invariant diffusions on invertible 
orientation scores (complex- valued functions on SE{2)) we put a lot of emphasis on the underlying differential ge- 
ometry in SE{2). All results straightforwardly generalize to the case of HARDI images, which can be considered 
as functions on R^ x S'^ embedded in SE{3). These rather technical results are summarized in Theorem [3] which 
answers all questions raised in the beginning of this section. Unfortunately, this theorem requires general differential 



^Or later in Subsectionjlsjto get some intuition in the choice of functions {ai}^^-^ and {D*-^ }^ 
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geometrical concepts such as principal fiber bundles, associated vector bundles, tangent bundles, frame-bundles and 
the Cartan-Ehresmann connection defined on them. These concepts are explained in full detail in 1631 (with a very 
nice overview on p. 386 ). 

The reader who is not familiar with these technicalities from differential geometry can skip the first part of the theorem 



while accepting the formula of the covariant derivatives given in Eq. (67), where the anti-symmetric Christoffel 



symbols are equal to minus the structure constants c^j = —Cj^ (recall Eq. (28)) of the Lie-algebra. Here we stress 
that we follow the Cartan viewpoint on differential geometry, where connections are expressed in moving coordinate 
frames (we use the frame of left-invariant vector fields {^i, . . . ^Aq} derived in Subsection [5] for this purpose) and 
thereby we have non- vanishing torsion]^ This is different from the Levi-Civita connection for differential geometry 
on Riemannian manifolds, which is much more common in image analysis. The Levi-Civita connection is the unique 
torsion free metric compatible connection on a Riemannian manifold and because of this vanishing torsion of the 
Levi-Civita connection V there is a 1-to-l relatior|^to the Christoffel symbols (required for covariant derivatives 
\/iV^ = diV^ + r^jdkV^ ) and the derivatives of the metric tensor. In the more general Cartan connection outlined 
below, however, one can have non- vanishing torsion and the Christoffels are not necessarily related to a metric tensor, 
nor need they be symmetric. 

Theorem 3 The Maurer-Cartan form uo on SE{3) is given by 

6 

^g{Xg) = ^{dA'\^,Xg)A,, XgeTg{SE{3)), (63) 

i=l 

where the dual vectors {dA^}^^i are given by (m) and Ai = It is a Cartan Ehresmann connection form on 

the principal fiber bundle P {SE{3),7r : SE{3) e = SE{3)/SE{3), SE{3)), where 7r{g) e, RgU ug, 
u^g e SE{3). Let Ad denote the adjoint action of SE (3) on its own Lie-algebra Te{SE (3)), i.e. Ad{g) = {Rg-iLg)^, 
i.e. the push-forward of conjugation. Then the adjoint representation of SE(3) on the vector space C{SE{3)) of left- 
invariant vector fields is given by 

Ad{g) =dnoAd{g)ouj. (64) 

This adjoint representation gives rise to the associated vector bundle SE{3) x ^ C{SE{3)). The corresponding 
connection form on this vector bundle is given by 

6 6 

cD = ^ ad{Aj) (g) d^^' = ^ c^j Au d^' (g) dA^ , (65) 

_ _ 6 

with ad = (Ad)^, i.e. ad{Aj) = ^^[A^^j] (X) dA\i50. p.265]. Then uj yields the following 6 x6-matrix valued 1 -form 

1=1 

Cj^{-) := -Cj{dA\ Aj) ^, j = 1, 2, . . . , 6. (66) 

on the frame bundle, ^63\ p.353,p.359], where the sections are moving frames R63\ p. 354]. Let {/ife}fc=i denote the 
sections in the tangent bundle E := {SE(3)^T{SE(3))) which coincide with the left-invariant vector fields {Ak}\^i' 
Then the matrix-valued 1-form given by Eq. ([66l) yields the Cartan connection given by the covariant derivatives 



Dxi^,,,iKim := DM7(i))(^l^(*)) 

k=l /c=l j=l (O/j 

k=l i,j,k=l 



^The torsion tensor Tv of a connection V is given by T^[X, Y] = VxY — VyX — [X, Y]. The torsion-tensor of a Levi-Civita 
connection vanishes, whereas the torsion-tensor of our Cartan connection V on SE{3) is given 

^In a Levi-Civita connection one has F^^ = FJ^ = | 9^^{9mk,i + 9mi,k — 9ki,m) with respect to a holonomic basis. 



31 



with a^{t) = 7*(^) ( ), for all tangent vectors -^l^(^) = 7^^) ^^l-yft) along a curve t ^ j{t) G 



-fit) 



i=l 



hit) 



SE{2) and all sections ii{^(t)) = Y f^k{l{t))- The Christoffel symbols in {67) are constant = — c^/., 

k=\ I I 

with c^^ the structure constants of Lie-algebra Te{SE{3)). Consequently, the connection D has constant curvature and 
constant torsion and the left-invariant evolution equations given in Eq. {42) can be rewritten in covariant derivatives 
(using short notation Vj := D^.): 



dtW{g,t) = E -a\W)A^W{g,t)+ E A^ {{D'' {W)){g,t)AjW){g,t) 

i—1 i^j — l 

= E -a\W)ViW{g,t)+ E Vi{{D'^{W))(g,t)VjW)(g,t) 

i—1 i^j — l 

W{g,0) = U{g) , for all ge SE{3),t> 0. 



(68) 



Both convection and diffusion in the left-invariant evolution equations ( [42] ) take place along the exponential curves 

6 

t E c'Ai 

lc,g{t) = g - e in SE{3) which are the covariantly constant curves (i.e. auto-parallels) with respect to the 

Cartan connection. In particular, ifa'^{W) = constant and ifD'^^{W) = (convection case) then the solutions are 



-t c'A, 

W{g,t) = U{g-e ). 



(69) 



The spatial projections Pr3 7 of these of the auto-parallel/exponential curves 7 are circular spirals with constant 
curvature and constant torsion. The curvature magnitude equals ||c^^^ ||~^ ||c^^^ x c^^^|| and the curvature vector 
equals 



(2) V ni^) 



^(2)1 



X ' X c 



(1) 



(70) 



where c = [c^ ^(? ^(? \ c^, c^, c^) = (iP^^ ; c^^^). The torsion vector equals T{t) = \ci • C2I n{t). 



Proof The proof is a straightforward generaUzation from our previous results fSO", Part II, Thm 3.8 and Thm 3.9] 
on the 6'£^(2)-case to the case SE{3). The formulas of the constant torsion and curvature of the spatial part of 
the auto-parallel curves (which are the exponential curv es) f ollow by the formula ( [tT] ) for (the spatial part x{s) of) 

Here we stress that s{t) = t ^/{c^^^^T~(^y^T~{^ 

x{s) and t{s) = 



9.1 



the exponential curves, which we will derive in Section 

is the arc-length of the spatial part of the exponential curve and where we recall that k,{s) 

(x(5) X x(s)). Note that both the formula d71| for the exponential curves in the next section and the formulas for 
torsion and curvature are simplifications of our earlier formulas | 39 , p. 175-177]. In the special case of only convection 



the solution Mj) follows by e^'^^^^^U{g) 
Ai = dn{Ai)7 



UetAUig), with A = - E-=i c'Ai and dn{A) = - E-=i c'Ai with 



9.1 The Exponential Curves and the Logarithmic IVIap expHcitly in Euler Angles 

The surjective exponent mapping is given by 

t EcM, _ J (cit,C2t,C3t , /)ifc^'^ =0, 
^'^^ ^ ^ \ {tC^^^ + ^-'p^'^ Qc^^^ + (tq-^ - ^H^) 1^2^(1) , / + -jj^Q + (l-cos(gt)) ^2) ^^^^ 

(71) 

where q = ||c^^^ || = ^/{c^yT{c^)^T{^ and Qx = c^^^ x x for all x G M^. 
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The logarithmic mapping on S'^(3) is given by: 



lo^ 




3 
i=l 


6 






3 


6 










Expressed in the first chart (using 


short notation 


= 4,7,/3,a)we have 





(72) 



q = arcsin y cos^ (^y^) sin^ /3 + cos"^ (^f ^ sin^(a + 7) ^^^^ 
c(^) := (c"^, c^, c^)^ = 2^1^ (sin /3(sina — sin 7) , sin /3(cosa + cos 7) , 2 cos^ sin(a + 7)) 

and 

c^i) = (c\ c^, c^)^ = X - ^ c(2) X X + r '(1 - (f ) cot(f )) c(2) X (c(2) X x). (74) 

Withc(l) = (C\C2,C3),C(2) = (C^C^C6). 

Throughout this article we will take the section a = (this is just a choice, we could have taken another section) in 
the partition x S'^ / ({0} x 5*0 (2)), which means that we will only consider the case 

(M) . p(2)\ _ (^1 ^2 3 4 5 6 N 

^ ' ~ ^ ~ v<-x,7,/3,a=05 ^x,7,^,05 <-x,7,^,05 ^x,7,/3,05 ^x,7,^,05 <-x,7,/3,0>'- 

Expressed in the second chart the section a = coincides with the section a = and along this section we again have 
( [74| ) but now with 

q = arcsin \Jcos^{j/2) sin^(;^) +cos^(/^/2) sin^(7) , ^^^^ 
c(2) = (c^,c^,c^)^ = ^T^^ (sin7 cos^(f ) , sin^^ cos^(|) , ^sin7 sin^^)^ , 

where we again used short notation c* •= ~ ^ Roughly speaking, c*^^^ is the spatial velocity of the exponential 
curve (fiber) and c*^^^ is the angular velocity of the exponential curve (fiber). 



10 Analysis of the Convolution Kernels of Scale Spaces on HARDI images 

It is notorious problem to find explicit formulas for the exact Green's functions : x S*^ of the left-invariant 



diffusions (44) on x 6'^. Explicit, tangible and exact formulas for heat-kernels on SE{3) do not seem to exist in 
literature. Nevertheless, there does exist a nice general theory overlapping the fields of functional analysis and group 
theory, see for example (681 (53, which at least provides Gaussian estimates for Green's functions of left-invariant 
diffusions on Lie groups, generated by subcoercive operators. In the remainder of this section we will employ this 
general theory to our special case where R^ x S'^ is embedded into SE{3) and we will derive new explicit and useful 



approximation formulas for these Green's functions. Within this section we use the second coordinate chart ( 15 ), as 



it is highly preferable over the more common Euler angle parametrization (12) because of the much more suitable 
singularity locations on the sphere. 

We shall first carry out the method of contraction. This method typically relates the group of positions and rotations 
to a (nilpotent) group positions and velocities and serves as an essential pre-requisite for our Gaussian estimates and 
approximation kernels later on. The reader who is not so much interested in the detailed analysis can skip this section 



and continue with the numerics explained in Chapter 12 
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10.1 Local Approximation of SE{3) by a Nilpotent Group via Contraction 



The group SE{3) is not nilpotent. This makes it hard to get tangible explicit formulae for the heat-kernels. Therefore 
we shall generalize our Heisenberg approximations of the Green's functions on SE{2), |33|, |65|,|4|, to the case 
SE{3). Again we will follow the general work by ter Elst and Robinson [ 68 1 on semigroups on Lie groups generated 
by weighted subcoercive operators. In their general work we consider a particular case by setting the Hilbert space 
]L2{SE{3)), the group SE{3) and the right-regular representation 1Z. Furthermore we consider the algebraic basis 
{^3, ^4, ^5} leading to the following filtration of the Lie algebra 

Qi = span{^3,^4, A} C 02 = span{^i,^2,^3,^4, A, A} = ^SE{3)) . (76) 

Now that we have this filtration we have to assign weights to the generators 

ws = W4 = = 1 and wi = W2 = wq = 2. (77) 

For example ws = 1 since A3 already occurs in Qi, wq = 2 since is within in ^2 and not in ^i. 

Now that we have these weights we define the following dilations on the Lie-algebra Te{SE{3)) (recall Ai = AilJ: 

Iqiic' A,) = 4 q^^ c' A,, for all d G M, 

=1 i=i 



and for < g < 1 we define the Lie product [A,B]q = ^~^[^q{A),^q{B)\. Now let {SE(3))q be the simply 
connected Lie group generated by the Lie algebra {Te{SE(3))^ ['i']q)- This Lie group is isomorphic to the matrix 
group with group product: 

^7^a) 'q (^'' = (X + • R^q^^q^^q2 ' Sq-i , R~^~ • i^y^,^, ) (78) 

where the diagonal 3 x 3-matrix is defined by Sq := diag{l, 1, q} and we used short-notation R;^^^ = Re^,jR^^ ^Re^,ar 
i.e. our elements of SO {3) are expressed in the second coordinate chart ( 15 ). Now the left-invariant vector fields on 
the group {SE{3))q are given by 

^i\g = ''Lg''lq)*^i^ i = 1,...,6. 

Straightforward (but intense) calculations yield (for each ^ = (x, R;^^^) G {SE{3))q ): 

A\\g = cos(^^q;) cos(^/3) dx + (003(7^') sin{aq^) + cos(a^^) sin0q) sin{^q)) dy + 

+q{s\n{aq^) sin(7g) — cos{aq^) cos{^q) sin(/3Q)) dz 
Al\g = — sm{aq'^) cos(/3^) dx + (cos(^^a) cos(7^) — sm{aq'^) sm{Pq) sin(7^)) dy-\- 

+^(sin(a^^) sm{Pq) cos(/3q) + cos(a^^) sin(7^)) dz 
Al\g = q~^ sm0q) dx - q~^ cos0q) sm{jq) dy + cos0q) cos{^q) dz 
Al\^ = -q~^ cos{aq^)tan0q) da + sin{aq) d^ + d^ 



: q ^ sin{aq^) tan(/3^) da + co^{aq^) d^ - "-^^ ^1 
da- 



6 

Now note that = 7g^M7g(^i). 7q(^j)] = lq^q'^''''^^[Ai, Aj] = g'^^^'^^^ "'^'^ 4 and thereby we have 

[A4,A2]q = q^A3, [A5,Ai]q = -q^As, [A5,A3]q = Ai, [Ae, Aijg = and [Ae, Asjg = -9^1. 

Analogously to the case q = I, {SE{3))q=i = SE{3) we have an isomorphism of the common Lie-algebra at the 
unity element Te{SE{3)) = Te{{SE{3))q) and left-invariant vector fields on the group {SE{3))q : 

{A, ^ A\ and A, ^ ^ \M, A,\q ^ \A\, A'^A . 
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It can be verified that the left-invariant vector fields satisfy the same commutation relations {19). 
Now let us consider the case q iO, then we get a nilpotent-group {SE{3))o with left-invariant vector fields 

^? = S,, ^0 ^ g^^ ^0 ^ ^g^ _ ^g^ ^ g^^ ^0 ^ _^g^ ^ g^^ ^0 ^ ^ . (80) 

10.1.1 The Heisenberg-approximation of the Time-integrated /c-step Contour Completion Kernel 

Recall that the generator of contour completion diffusion equals A3 + 1)^^((^4)^ + (^5)^). So let us replace the true 
left-invariant vector fields {^^jf^g on SE{3) = {SE{3))q=i by their Heisenberg-approximations {A^}^^^ that are 

given by ( 80 ) and compute the Green's function Pt^~ ' ~ on {SE{3))o (i.e. the convolution kernel which yields 



the solutions of contour completion on {SE{3))o by group convolution on {SE{3))o). For < D^"^ « 1 this kernel 
is a local approximation of the true contour completion kernej^p^^"^'^ ~^ , onR^ y\ S"^: 



12(a3-(l/2)z/3)2+z2;32 12 (y + ( 1 /2) z^) ^ 7^ 

where n(/3, 7) = R^^^^yR^^ ^e^ = (sin /3, — sin 7 cos /3, cos 7 cos P)^ . The corresponding /c-step resolvent kernel on 
the group {SE{3))o is now directly obtained by conditional integration over tim^ 



kk-5 . 12(x-(l/2)z/3)2+z2/32 12 (y + ( 1 /2) z^) ^ 7^ 

ifz < Oand(x,^) / (0,0). 



10.1.2 Approximations of the Contour Enhancement Kernel 

Recall that the generator of contour completion diffusion equals D^^{A3)'^ + I)^^((^4)^ + (^5)^). So let us replace 
the true left-invariant vector fields {Ai}^^^ on SE{3) = {SE{3))q^i by their Heisenberg-approximations {A^}^^^ 



given by ( 80 ) and consider the Green's function Pt ' ~ on {SE{3))o'. 



Now since the Heisenberg approximation kernel ' ' ^ ^ is for reasonable parameter settings (that is < 

44 _D^'^ D"^^ ■ (SE(2)) 

^33 << 1) close to the exact kernel p^ ' ' ^ ^ we heuristically propose for these reasonable parameter settings 
the same direct-product approximation for the exact contour-enhancement kernels on x : 

Pt {x,y,z,n{l3,-f)) ^ N{D ,D ,t)p^ ^ '\z/2,x, -p^ ' (^/2, -2/, 7) , (83) 

where 

N{D^^, D^^, t) ^ -^V^tVw^VD^^ (84) 
8v 2 



^The su pers cript for the kernel is actually so in the superscript-labels, for the sake of simplicity, we only mention the non-zero coefficients 
D'^,ai of (|44|. 

^Note that the delta distribution S{s — z) allowed us to replace all s by 2; in the remaining factor in (81} which makes it easy to apply the 
mtegration = p^^ ' r(f;fc,A)dt. 
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takes care of Li (R^ x S*^) -normalization and with 



— JJ ,JJ / p\ 1 

Pt [^jI/t^) = 47r^2i:)44^33 ^ 



£)33 £,44 



r-,33 7-^44 



47rt2D44^33C 



1 If 9'^ e^jy-i-x sine + y cose))^ \^ l I ( (x cos g + y sin 0) - x) I - 

^2YVo44^ 4(l-cos(a))2D33 ^ ^o44i^33| 2(l-cos0) | 



, if ^7^0, 



1 p 4tc2 

^ 47rt2D44D33^ 



X- ) I ly|- 

033 ; ^i544o33 _ Q 



(85) 

which are reasonably sharp estimates of hypoelliptic diffusion on SE{2), with ^ < c < v^, for details see (301 ch 



5.4]. For the purpose of numerical computation, we simplify ' (x, 6) in ( 85 ) to 



Vt 



47rt2i:)44^33 



fey , g/2 \2X 2 

^2 , V 2 +tan(g/2) ] 1 /^-ccg , 6/2 

^44+ ^33 +^^44^33 ^ 



2 "^tan(6»/2) 



where one can use the estimate 



tan%%) ^ i-(1v2l) 1^1 <fo^^ a^oi^ numerical errors. 



10.2 Gaussian Estimates for the Heat-kernels on SE{3) 

In (261 ch:6.2] it is shown that the constants Cs, C4 are very close and that a reasonably sharp approximation and 
upperbound of the horizontal diffusion kernel on x /S^ is given by 

'°>(y,n(^,7)) « (,,,,,^6-^^^^^^ (86) 

with weighted modulus 



liy,^nj|D33,D44 .- y ^33^44 + + 1^^D33" ^ ^44 J ^^/^^ 

where n = n(/3, 7) G 5*^ and where we again use short notation := c^^i(x, ^ q), /c = 1, . . . , 6. Recall from 



Section 9.1 that these constants are computed by the logarithm (72) on SE{3) or more explicitly by (75 ) 



10.3 Analytic estimates for the Green's functions of the Cost Processes on >^ S'^ 



In Appendix \B\ we have derived the following analytic approximation for the Green's function k 



of the 



Hamilton- Jacobi equation (136) and corresponding cost-process explained in Section 8.3 Eq. (62) on R^^ x S 



IT] 



|P(y,a = 0,/3,7)P 



^2 



(88) 



for sufficiently small time t > 0, with weights Wi = W2 = w/^ = = 1, w^, = wq = 2 and 2 > C > and where 
the (c^, . . . , c^) are given by y5\ and where we have set n = n{p^ 7) and := c^^i(y, R~ ^ q). Here we recall 
that D^^ tunes the spatial erosion orthogonal to fibers, D^^ tunes angular erosion and r] G (|, 1] relates homogeneous 
erosion r] = \io standard quadratic erosion r] = 1. 



Now analogously to our approach on SE{2) in (30l Ch:5.4] we use the estimate 

\a\ + \b\>^W + W>^{\a\ + \b\) 
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a, 6 G M to obtain differentiable analytic local approximations: 



( \c'\' + V 



^2|2 |^4|2 I |^5|2\2 |^3|2 \ 2?7-l 

\ ^44 ) -r ^TT^ ] ,7/^(^2'-'-] 



(89) 



and for 7^ = ^ we obtain the flat analytic local approximation (that arises by taking the limit r] 



(y,n) ?^ { Y V ^ / D^^D^^ - 



(90) 



else 



See Figure 13 for glyph visualizations (recall Definition [T]) of the erosion and corresponding diffusion kernel 

on the contact manifold (R^ >^ 6*^, d^^). In order to study the accuracy of this approximation formula for the case of 





Figure 13: Left: Diffusion kernel in Eq. (91 ), where angular diffusion takes place over the full S'^ sphere and spatial 
diffusion simultaneously takes place in the plane G ' ' -orthogonal to As. Right: Erosion kernel in ( 89 ). 



angular erosion only (i.e. D^^ = 0) we have analytically computed 



m(^,7) 



I3ll=0,I?44 = i 2 



_ 1 (a^((c^)^+?j^))^+cos-^ ^jd^iic^fHc'f))' 

~ 4 (c4)2 + (c5)2 



dtk^ 



where we used the following identities 



(92) 



\A^U\' + \A,U\' = \djU\' + (cos^)-2|S?/7p 

\AiU\^ + \A2U\^ = |S^/7p + |S^/7p + \dzU\^ - lAsU]^ for all U e C\R^ x S^). 
Ideally m = 1 since then the approximation is exact. For relevant parameter settings we indeed have m ^ 1 up to 



5 -percent ^oo -errors as can be seen in Figure 14 



11 Pseudo Linear Scale Spaces on x 

So far we have considered anisotropic diffusions aligned with fibers and erosions orthogonal to the fibers. As these 
two types of left-invariant evolutions are supposed to be alternated, 
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44 



Figure 14: This figure shows that our approximation (89) is rather accurate for relevant parameter settings /3,7 G 
[~f 5 f ]• have computed the surface m = m{l3^j) where m(;^, 7) is given by (92). The approximation is exact iff 
m = 1, which is the case along the lines ^ = and 7 = 0. Relative errors are smaller than 0.05 percent. 



where g'^^ are the components of the inverse metric, which in case of a diagonal metric tensor simply reads = 
g^-^ = D^^, Hke in (52) where D^^ D^^ and where denotes the viscosity solution operator U ^ 



for the erosion Hamilton- Jacobi equation ( |52| ). the natural question arises is there a single evolution process that 
combines erosion/dilation and diffusion. Moreover, a from a practical point of view quite satisfactory alternative to 
visually sharpen distributions on positions and orientations is to apply monotonic greyvalue transformations (instead 
of an erosion) such as for example the power operator Xp 

fe(/7))(y,n) = (C/(y,n))P, 

where p > I, where we also recall the drawback illustrated in Figure |4] Conjugation of the diffusion operator with a 
monotonically increasing grey-value transformation x • 



(93) 



is related to simultaneous erosion and diffusion. For a specific choice of grey- value transformations this is indeed the 
case we will show next, where we extend the theory for pseudo linear scale space representations of greyscale images 
(3611 to DW-MRI (HARD! and DTI). 

Next we derive the operator o Q^'^CV) ^ ^ more explicitly using the chain-law for differentiation: 

Aj{x{V{;;t))) =x'{V{;;t))MV{;;t)), 
dtX{V{;;t)) =x'{V{;;t))dt{V{;;t)), 
AiAjixiVi; ; t))) = X"{V{; ■,t))Ai{V{; -,1)) Aj{V{; t)) + x'{V{; ■,t))AAj{V{; ; t)). 

Consequently, if we set = x{V{-,-,t)) ^ V{-,-,t) = x~HW{-,-,t)) and 



Wi;;t) 



D=diag{Dll,D",D33,D",D",0} 



,a=0 



Wi;;0), 



see Figure 15 then V satisfies 



5 5 

dtV{y,n,t) = Y,^D" {{A\^/V) (y, n, + ^^y, E ^" (A|(y,„) ny,n,t))' 
V{y,n,0)=x-HU{y,n)). 



(94) 



where fi{V{y, n, t)) = \i(x{v{y'n't))) ■ i^i^iy, n, t j) = C constant we achieved that (|93| coincides with 

a simultaneous erosion/dilation and diffusion with 
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X 



X 



vo,-,o) 



clilation<S^ 
diffusion 



Figure 15: The commutative diagram of simultaneous left-invariant dilation and diffusion along the fibers (D'^'^ 



D^^) given by Eq. (94) with fl{V{y^ n, t)) = C illustrates that conjugation of specific grey- value transformations (96) 



with left-invariant diffusion is equivalent. 



Figure 16: Left: Graphs of the grey- value transformations, Eq. (96), that are solutions xc • [0, 1] [0, 1] of Eq. (95), 
forC = —8,— 4,— 2,— 1,0, 1,2, 4, 8, 1361 (depicted from left to right). The concave solutions correspond to diffusion 



and erosion whereas the convex solutions yield convection and diffusion. Right: Graphs of the convex xc of Eq. ( [95 
depicted together with their inverse Xc^^ Eq (97 ) , for C = —8, —4, —2, —1 



This means we have to solve the following ODE-system 



x"{i)-Cx'{i) = ^ /e[o,i] 

X(0) = 0andx(l) = l, 



where / G [0, 1] stands for intensity where in particular we set / := 
The unique solutions of ([95]) are given by 



y(y,n,0)-min(y,„) y(y,n,0) 
maxy,n y(y,n,0)-min(y^n) V'(y,n,0) 



(y,n) G 



Xc{I) = 



if C 7^0, 
ifC = 0, 



(95) 
(96) 



so that V{y, n, t) is the solution of an evolution (44) where the generator is a weighted sum of a diffusion and ero- 
sion/dilation operator. The inverse of % is given by 



iln(l + (e^-l)/) if C 7^0, 
/ ifC = 0. 



(97) 



The drawback of this intriguing correspondence, is that our diffusions primarily take place along the fibers, whereas 
our erosions take place orthogonal to the fibers. Therefore, at this point the correspondence between pseudolinear 
scale spaces and hypo-elliptic diffusion conjugated with ^ is primarily useful for simultaneous dilation and diffusion 
along the fiber^ i.e. the case 

D = diag{0, 0, D^\D^^, D^^, 0} and C> 0. 



In general one does not want to erode and diffuse in the same direction. 
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Remark 4 Finally, for C ^ 2 the operator in Eq. {93) is close to the operators applied in Figure^and Figure^on 
DTI and HARDI data of the brain. The difference though is that ([9J]) conjugates a diffusion operator with greyvalue 
transformations, whereas in Figure^and^we applied Xp o e^^^^"^^^ o Xp "^ith p = 2. 



12 Implementation of the Left-Invariant Derivatives and x -Evolutions 



In our implementations we do not use the two charts (among which the Euler-angles parametrization) of 5*^ because 
this would involve cumbersome and expensive bookkeeping of mapping the coordinates from one chart to the other 
(which becomes necessary each time the singularities ([14]) and ([17]) are reached). Instead we recall that the left- 
invariant vector fields on HARDI-orientation scores U : SE{3) M, which by definition (recall Definition [s]) 
automatically satisfy 

U{RR,^J = U{R), (98) 
are constructed by the derivative of the right-regular representation 

AU{g) = {dn{A,)U){g) = lim ^ = lim ^---^ ^ , 

h^O fl hlO Zfl 

where in the numerics we can take finite step-sizes in the righthand side. Now in order to avoid a redundant computa- 
tion we can also avoid taking the de-tour via HARDI-orientation scores and actually work with the left-invariant vector 
fields on the HARDI data itself. To this end we need the consistent right-action 91 of SE{3) acting on the space of 
HARDI images L2(M^ x 5'^). To construct this consistent right-action we formally define S : L2(M^ x S'^) ^ H, 
where H denotes the space of HARDI-orientation scores, that equals the space of quadratic integrable functions on 
the group SE{3) which are a right-invariant, i.e. satisfying ([98]) by 



{SU){x,R) = U{x,R) = U{x,Re,). 

This mapping is injective and its left-inverse is given by {S~^U){x^ n) = U (x, R^), where again R^ G 5*0 (3) is some 
rotation such that i^nC^ = n. Now the consistent right-action 9^ : SE{3) 5(L2 (M^ x S'^)), where 5(L2(M^ x 5^)) 
stands for all bounded linear operators on the space of HARDI images, is (almost everywhere) given by 

(%,i?)^)(y,n) = (^-^ o7^(,,^) o^[/)(y,n) = U{RnX ^ R^Re,). 

This yields the left-invariant vector fields (directly) on sufficiently smooth HARDI images: 

AMy,n) = imA.)U)iy,n) = (^e^^,/7)(y n) - ;7(y,n) ^ ^.^ (91...,/7)(y,n) - (fn,-...^)(y,n)_ 

h^o h h^o 2h 

Now in our algorithms we take finite step-sizes and elementary computations (using the exponent ( jTT] )) yield the 
following simple expressions for the discrete left-invariant vector fields, for respectively central. 



AlU(y, n) ^ AlU{y, - U(y+hR^e^ ,n)-U(y-hR^e^ ,n) 

A2U{y,n)^AW{y,n 



U(y+h Rney , n)-U(y-h R^ey , n) A4U [y , U) ~ A4U [y, U) .- ^ 

2h ' A TT/ \ \CTT( \ U{y , Rn R^^^h ^z)-U{y , Rn R^^^-h ^z) 
U(y+hRne^ ,n)-U(y-hRne^ , n) A5U (y, II) ^ A^U^y, II) := ^ ^ ^ 

(99) 



2h 



^3f/(y,n)«^i[/(y,n) 
forward, 

A^Uiy, n) « A{Uiy, n) := ^(y+^'^ne ")-r/(y,n) ujy H,^„,..)-uiy 

u{y+hR^ey ,n)-u(y,n) Ji4U [y,n) {y,n) .- 



A2U{y,n)^AlU{y,n) 
A3U{y,n)^AiU{y,n) 



and backward. 



h 

AiU(y,n)^A\U{y,n) •= U(y ,n)-U(y-h R^e^ ,n) 

A2U{y,n)^AlU{y,n) 
A3U{y,n)^A'sU{y,n) 



U{y+hRneJ',n)-U{y,n) ' A t/(y, n) ^ A^U (y, II) 



U{y,n)-U{y-hRney,n) ' ^4^7 (y, n) ^ A\U (y, II) 
^(y,n)-^(y-^Hne. ,n) ' A->U (j , II) ^ AIU (y, II) 









R,y!h e,)-t/(y,n) 




h 


U{y,n)- 


U{y,Rn R,^,-h ez) 






U(y,n)- 


U(yXn R,y,-h ez) 



(100) 



(101) 



h 
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left-invariant finite differences. The left-invariant vector fields {^1,^2,-4.4,^5} clearly depend on the choice of 
Rn G 5*0 (3) which maps Rn^z = n. Now functions in the space H are a-right invariant, so thereby we may assume 
that R can be written sls R = R^^^^R^y^is, now if we choose Rn again such that Rn{i3,-f) = ^e^,7^e^,/3^e^,Q!=ao=o 
then we take consistent sections in SO{3)/SO{2) and we get full invertibility o S = S o = X. 

Our evolution schemes, however, the choice of representant Rn is irrelevant, because they are well-defined on the 
quotient R3 x 5^ = SE{3)/{{0} x S0{2)). 



In the computation of (99) one would have liked to work with discrete subgroups of SO (3) acting on S'^ in order to 
avoid interpolations, but unfortunately the platonic solid with the largest amount of vertices (only 20) is the dodeca- 
hedron and the platonic solid with the largest amount of faces (again only 20) is the icosahedron. Nevertheless, we 
would like to sample the 2 -sphere such that the distance between sampling points should be as equal as possible and 
simultaneously the area around each sample point should be as equal as possible. Therefore we follow the common 
approach by regular triangulations (i.e. each triangle is regularly divided into (o + 1)^ triangles) of the icosahedron, 
followed by a projection on the sphere. This leads to No = 2 + 10(o+ 1)^ vertices. We typically considered o = 1, 2, 3, 
for further motivation regarding uniform spherical sampling, see | 39, ch.7.8.1]. 



For the required interpolations to compute ( [99] ) within our spherical sampling there are two simple options. Either one 
uses a triangular interpolation of using the three closest sampling points, or one uses a discrete spherical harmonic 
interpolation. The disadvantage of the first and simplest approach is that it introduces additional blurring, whereas the 
second approach can lead to overshoots and undershoots. In both cases it is for computational efficiency advisable to 
pre-compute the interpolation matrix, cf. |21 , ch:2.1], ||4TJ p. 193]. See Appendix|F| 



12.1 Finite difference schemes for diffusion and pseudo-linear scale spaces on >^ S'^ 

The linear diffusion system on x /S^ can be rewritten as 

dtWiy,n,t) = (D"((A)2 + (A2f)+D^^iAs)^+D^^As2) Wiy,n,t) 



l^(y,n,0) = C/(y,n), 



(102) 



This system is the Fokker-Planck equation of horizontal Brownian motion on x if D^^ = 0. Spatially, we 
take second order centered finite differences for (^1)^, (^2)^ and (^3)^, i.e. we applied the discrete operators in the 



righthand side of (99 ) twice (where we replaced 2h hto ensure direct-neighbors interaction), so we have 



iiAsfWKy,n,t) « Wjy + hR.e.,n,t) - 2Wiy n,t) ^Wjy - hR.e.,n,t) ^^^^^ 

For each G N U {0}, we define the vector = {w^ f.)k=i...Nojex ^ R^^^°, where k enumerates the number of 
samples on the sphere, where y enumerates the samples on the discrete spatial grid X = {1,...,A^} x {1,...,A/'} x 
{1, . . . , A^} and where q enumerates the discrete time frames. Rewrite ( 103 ) and ( |102| ) in vector form using Euler- 
forward first order approximation in time: 



W^+l = (/-At(J5,+JM3))w^ 

the 

the spatial increments block-matrix. 



where J52 G ^^^^oxn^n^ (denotes the angular increments block-matrix and where Jj^s GG ^^^^oxn^No (denotes 



12.1.1 Angular increments block-matrix 

The angular increments block-matrix equals in the basic (more practical) approach, cf. lSTIl . of (tri-)linear spherical 
interpolation 

Js,=h-^D^\B^^B,), (104) 
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with angular stepsize ha > and with B4 and are matrices with 2 on the diagonal and where for each column the 
of off-diagonal elements is also equal to 2. See Appendix |F| Eq. ( |164| ), for details. 



In the more complicated discrete spherical harmonic transform interpolation approach, cf. |[34l Ch:7] the angular 
increments block-matrix is given by 

(105) 

with regularization parameter treg > and where Nsh represents the number of spherical harmonics and with 
nsH ^ TVo -matrix 

M = [Ml] = l^Y^gink)], with^i) = LVT^J and m(j) = j - {l{j) f - - 1 

with C = Yjjli \^m(j) 0) P where the diagonal matrix A = dmg{5s-2 (ni), . . . , (5^2 (hat^)} contains discrete 
surface measures Ss^ (n/e) (for spherical sampling by means of higher order tessellations of the icosahedron) given by 

6s2{nk) = ^ A{ni,nj,nk), (106) 

where i ^ j means that and Uj are part of a locally smallest triangle in the tessellation and where the surface 
measure of the spherical projection of such a triangle is given by 



A{ni,nj,nk) = 4 arctan(Y^tan(si_^-/e/2) tan((si_^-/e - Sij)/2)t8in{{sijk - Sik)/2)tdin{{sijk - Sjk)/2)) 
with Sijk = ^{sij + Sik + Sjk) and Sij = arccos(ni • n^). 



12.1.2 Spatial increments matrix 

In both approaches the spatial increments block matrix equals Jj^s given by Eq. ( |164| ) in Appendix |F| 



12.1.3 Stability bounds on the step-size 

We have guaranteed stability iff 

\\I - At{Js2 ^ Jms)\\ =sup||^||^i||(/- At(J52 + Jm3))w|| < 1, 
which is by / - At( J52 +75^3) = {ul - At Jrs) + ((1 - tj)I - At J52), with u e (0, 1), the case if 

||/ - ^ Jr3 II < 1 and ||/ - ^ J52 II < 1 ^ II (^/ - At Jms) + ((1 - - At J52) || < + (1 - z.) = 1. (107) 



Sufficient (and sharp) conditions for the first inequality are obtained by means of the Gerschgorin Theorem {43l : 

^ 4^11^2^33 - (108) 

Application of the Gerschgorin Theorem li43J on the second equality gives after optimization over the following 
global guaranteed- stability bound 

< 4^11^2^33 I (109) 
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in the linear interpolation case and 
At < 



4Dll+2D33+D44/j,2_ 



At < 



4Dii+2D33+D44/j^2^ 



if treg - L{L ^ 1) < 1 
if treg ■ L{L ^ 1) > 1 



(110) 



in the spherical harmonic interpolation case. For details, see respectively | 21 1 and 1341 . 



The pseudo-linear scale spaces are implemented using the above together with Eq. (93 ), and Eq. (96 ), Eq. (97 ). 

12.2 Finite difference schemes for Hamilton- Jacobi equations on >^ S'^ 



Similar to the previous section we use a Euler-forward first order time integration scheme, but now we use a so-called 
upwind-scheme, where the sign of the central difference determines whether a forward or central difference is taken. 



At 



a-'^ (y, n) U(y ,n)-U(y-h R.e^ ,n) ^ ^+,1 ^ 



iy(y,n,t + At) = Ty(y,n)+ 

+ (a-'^ (y, n) ' ""^-""^-^ ^"^^ ' + a+'^ (y, n) "^"^^^ ' ) ' | + 



Ujy+hRnea: ,n)-L7(y,n) 
h 



+ 



where the functions a are given by 



'"^(y, n) = max < 0. 



U (y+ hRn^x ,n)-U (y- h, HnCx , n) 



a-'2(y,n) = maxjo, ^(y+^ ^ne, , nKL7(y-^ i^ne, , n) 



a+'^(y,n) = minjo, 
a^' (y, n) = mm < 0, ■ 



2h j ' 

} 



+''(y,n) 


= min < 




+''(y,n) 


= min < 




-•"(y,n) 


= max 


r 

0, 


-•'(y,n) 


= max 





2h 

U(y,Rn Rex,h e^)-U(y,Rn R^^ 



2h 



U(y,Rn R,y,h e,)-U(y,Rn R,y,-h e.) ^ 
2h J • 



12.3 Convolution implementations 



The algorithm for solving Hamilton- Jacobi Equations Eq. 136 ( 52 ) by x S'^ dilation/erosion with the corresponding 
analytic Green's function (89) boils down to the same algorithm as x 6*^ -convolutions with analytic Green's 



functions for diffusion. Basically, the algebra (+, •) is to be replaced by respectively the (max, +) and (min, +)- 
algebra. So the results on fast efficient computation, using lookup-tables and/or parallelization cf. 1601 , apply also to 
the morphological convolutions. The difference though is that hypo-elliptic diffusion kernels are mainly supported 
near the z-axis (thereby within the R^ x S'^ convolution after translation and rotation near the A3 -axis), whereas the 
erosion kernels are mainly supported near the xy-plane, but the efficiency principles are easily carried over to the 
morphological convolutions. 

Next we address two minor issues that arise when implementing the convolution algorithms |34, ch:8.2] and I6QII : 



1 . In the implementation with analytical kernels expressed in the second chart, such as ( 83 ) one needs to extract the 
second chart Euler angles {^^ 7) from each normal n in the spherical sampling, i.e. one must solve n{/3^ 7) = n. 
The solutions are 



(^,7) = (sign(ni) arccos(sign(n2^n| + n|)), -sign(n3) arcsin ( ^ j) (111) 
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with sign(x) = 1 if x > and zero else, and with sign(x) = 1 if x > and zero else. 

2. X 5^ -convolution requires computation of R^y , with v = y — Theorem[l] For the sake of computation 
speed this can be done without goniometric formulas: 

1 / {'ri'^)'^v^ + n^in?{n^ — l)v'^ + n^{n^n^v^ — \/ {n^Y + {n^Y^l — (n^yv^) \ 

) ^^"^ ) \ ^l-(n3)VK)' + (n2)2(ni^i + n^^^) ^ ^3^3((^i)2 ^ (^2)2) ^ 

if (n^, n^) 7^ (0, 0), where i?n is indeed a rotation given by 

W = 2A2 (n^)'^ ' + n^n^(n ^ - 1)^^ + n^n 'n^v^ + VK)^ + (n^) Vl - (n^) V) 

if (n\n^) ^ (0,0) that maps (0,0,1)^ onto (n\n^,n^) G 5'^. If (n^n^) = (0,0) then we set i?n=(o,o,n3) = 
sign(n^) /. 



13 Adaptive, Left-Invariant Diffusions on HARDI images 



13.1 Scalar Valued Adaptive Conductivity 



In order to avoid mutual influence of anisotropic regions (areas with fibers) and isotropic regions (ventricles) one can 
replace the constant diffusivity/conductivity D^^ by 



D^^A^A^ ^ D^^As o e o As 



(112) 



in the generator of the left-invariant diffusion (44 ) (where a = and D = diag{0, 0, D^^^ 0, 0}). This now yields the 
following hypoelliptic diffusion system 



W^(y,n,0) =/7(y,n), 



(e ^3^)(y,n,t), 



(113) 

Here could also choose to adapt the diffusivity by the original DW-MRI data : x S' ^ at time 0, so that the 

diffusion equation itself is linear, whereas the mapping U ^ ^t{U) := VK(-, t) in Eq. (113) included is well-posed 
and nonli near In the experiments, however we extended the standard approach by PeronaMalik [58 1 to x S'^ and 
we used (113) where both the diffusion equation and the mapping U ^ W{-^t) are nonlinear. The idea is simple: 
the replacement sets a soft- threshold on the diffusion in ^3 -direction, at fiber locations one expects |^3/7(y, n)| to 
be small, whereas in the transition areas between ventricles and white matter, where one needs to block the fiber 
propagation by hypo-elliptic diffusion, one expects a large \AsU{y^n)\. For further details see t he sec ond author's 
master thesis | 21 1. Regarding discretization of (|113|) in the finite difference schemes of Subsection 1 1 2. 1 1 we propose 



^3Ty(-,t)))(y,n) 



D^^ (y+ |h, n, t) • AW{y+ |h, n, t) - D^^ (y- |h, n, t) • A'sW{y- |h, n, t) 



h 



where h is the spatial stepsize and where h = hRnez for notational convenience and where 



>c{|^|w(y,n,t)|,|^|w(y,n,t)|} 



and where we recall Eq. ([99]), Eq. ( |1Q0| ) and Eq. ( |101| ) for central, forward and backward finite differences. 
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13.2 Tensor- valued adaptive conductivity 



For details see (34] and f26^, Ch:9], where we generalized the results in |[3T]|40l. Current implementations do not 
produce the expected results. Here the following problems arise. Firstly, the logarithm and exponential curves in 
Section 9.1 are not well-defined on the quotient x S'^ and we need to take an appropriate section through the 
partition of cosets in x 5^. Secondly, the. 6x6 Hessian in (261 Eq.l03,Ch:9] has a nilspace in ^e-direction. 



14 Experiments and Evaluation 



First we provide a chronological evaluation of the experiments depicted in various Figures so far. 



In Figure Hand Figure [3] one can find experiments of R^ x 5^ -convolution implementations of hypo-elliptic diffusion 
("contour enhancement", cf. Figure 10), Eq. (49), using the analytic approximation of the Green's function Eq. ( [83] ). 
Typically, these R^ x 6*^ -convolution kernel operators propagate fibers slightly better at crossings, as they do not suffer 
from numerical blur (spherical interpolation) in finite difference schemes explained in Section 12.1 Compare Figure[2] 
to Figure |5] where the same type of fiber crossings of the corona radiata and corpus callosum occur. Here we note 
that all methods in Figure [5]have been evaluated by finite difference schemes, based on linear spherical interpolation, 
with sufficiently small time steps At, recall Eq. ( |109| ) using 162 spherical samples from a higher order tessellation of 
the icosahedron. In Figure [5]the visually most appealing results are obtained by a single concatenation of diffusion 
and erosion. Then in Figure]4] we have experimentally illustrated the short-coming of sharpening the data by means 
of squaring in comparison to the erosions, Eq. (52). Experiments of the latter evolution, implemented by the upwind 
finite difference method described in Section [T^ has been illustrated for various parameter settings in Figure [8] (also 
in Figure [5]). The possibilities of a further modification, where erosion is locally adapted to the data, can be found in 
Figure |9] Then in Figure [T0| and Figure [TT] we have depicted the Green's functions of the (iterated) direction process 
(hypo-elliptic convection and diffusion) that we applied in a fiber completion experiment in Figure [12 



The experiments in Figure [5] (implemented in Mathematica 7) take place on volumetric DTI-data-sets 104 x 104 x 10 
where we applied ([2]) using 162 samples on the sphere (tessellation icosahedron). In Figure 17 we illustrate the 
enhancement and eroding effect in 31), of the same experiments as in Figure [5] where we took two different viewpoints 
on a larger 3D-part of the total data-set. 



In the remaining subsections we will consider some further experiments on respectively diffusion, erosion and pseudo- 
linear scale spaces. 



14.1 Further experiments diffusions 



In Figure 18 we have applied both linear and nonlinear diffusion on an noisy artificial data set which both contains an 
anisotropic part (modeling neural fibers) and containing a large isotropic part (modeling ventricles, compare to Fig. [5]). 
When applying linear diffusion on the data-set the isotropic part propagates and destroys the fiber structure in the 

anisotropic part, when applying nonlinear diffusion with the same parameter settings, but now with D = e 
and suitable choice of K > one turns of the propagation (diffusion-flow) through the boundary between isotropic and 
non-isotropic areas where (y, n |^3(y, n)p is relatively high. Consequently, both the fiber-part and the ventricle 
part are both appropriately diffused without too much influence on each other, yielding both smooth/aligned fibers and 
smooth isotropic glyphs. 
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Figure 17: Same experiments as in Figure [5] but now visualized with several slices in 3D. First row: Input DTI- 
data set. Second row: Output squared linear diffusion on squared data-set. Third row: Output erosion applied to the 
diffused dataset depicted in the second row. For parameter settings, see Figure [5] 



Input U = Ui + U2 




Output ^t(U) linear diffusion Output nonlinear diffusion 




Figure 18: Adaptive diffusivity/conductivity based on the data, Section 13.1 Top row: Artificial 15 x 15 x 15 x 162- 
input data that is a sum of a noisy fiber part and a noisy isotropic part. For the sake of visualization we depicted these 
parts separately. Bottom row: Output of just linear diffusion (without grey value transformations, nor Eq. (|4])), t = 1, 
D^^ = 1, D^^ = 0.04. Output nonlinear diffusion with t = 1, D^^ = 1, D^^ = 0.015, K = 0.05, At = 0.01. 
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input 



7] = 0.5 7] = 0.6 7] = 0.7 7] = 0.8 7] = 0.9 r}=1.0 




Figure 19: The effect onr] e [^,1] on angular erosion Eq. (52), D^^ = 0.4, D^^ = and t = 0.4, At = 0.02. Left: 
original glyph, right eroded glyphs for = 0.5, . . . , 1.0. 





Figure 20: Morphological kernel implementation of the Hamilton- Jacobi equations Eq. (52) by x -erosion, 
Eq. ([39]), 7/ = 1, D^^ = 0, D^^ = 0.4, t = 1, Left: input and output for No = 162, Right: same for No = 362. 



14.2 Further experiments erosions 



Erosions yield a geometrical alternative to more ad-hoc data- sharpening by squaring, see Fig}4]for a basic illustration. 
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where we fixed D"^"^ 0.4, D^^ and 



We compared angular erosions for several values of 7^ G [0, 1], in Figure 

t = 0.4 in Eq. ([52]). Typically, for r] close to a half the structure elements are more flat, Eq. ( [89] ) and Eq. ( [9Q| ) leading 
to a more radical erosion effect if time is fixed. Note that this difference is also due to the fact that D^^ in Eq. ( 52 ) has 
physical dimension [^"^^"^^j^jy and D""^ has physical dimension l^^i. 



[Time] 



We have implemented two algorithms for erosion, by morphological convolution in R'^ x S erosions Eq. ( 39 ) with 



analytical kernel approximations Eq. ( |89| ) and finite difference upwind schemes with stepsize t > 0. A drawback of 
the morphological convolutions is that it typically requires high sampling on the sphere, as can be seen in Figure [20 
whereas the finite difference schemes also work out well for low sampling rates (in the experiments on real medical 
data sets in [5] we used a 162-point tessellation of the icosahedron). Nevertheless, for high sampling on the sphere the 
analytical kernel implementations were close to the numerical finite difference schemes and apparently our analytical 
approximation Eq. (90) approximates the propagation of balls in y\ S'^ close enough for reasonable parameter- 
settings in practice. Upwind finite difference schemes are better suited for low sampling rates on the sphere, but for 
such sampling rates they do involve inevitable numerical blur due to spherical grid interpolations. Typically, erosions 
should be applied to strongly smoothed data as it enhances local extrema. In Figure [5j lower left corner one can see the 
result on a DTLdata set containing fibers of the corpus callosum and corona radiata, where the effect of both spatial 
erosion (glyphs in the middle of a "fiber bundle" are larger than the glyphs at the boundary) and angular erosion glyphs 
are much sharper so that it is easier to visually track the smoothly varying fibers in the output dataset. Typically, as can 
be seen in Fig.[5j applying first diffusion and then erosion is stable with respect to the parameter settings and produces 
visually more appealing sharper results then simultaneous diffusion and dilation in pseudo-linear scale spaces. 



14.3 Further experiments pseudo-linear scale spaces 



Pseudo-linear scale spaces, Eq [94] combine diffusion and dilation along fibers in a single evolution. The drawback is 

we 



a relatively sensitive parameter C > that balances infinitesimally between dilation and diffusion. In Figure [21 
applied some experiments where we kept parameters fixed except for C > 0. The generator of a pseudo-linear scale 
space consists of a diffusion generator plus C times an erosion generator and indeed for C = 4 we see more elongated 
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Figure 21: Pseudo-linear scale space on x 5^ applied on the initial DTI-dataset in Figure [s] t = D^^ = 1, 
D^^ = 0.04, At = 0.01, top: C = 2, bottom: C = 4. 



glyphs and more apparent crossings than for the case C = 2. 



15 Conclusion 



For the purpose of tractography detection and visualization of biological fibers. Diffusion- Weighted MRI-images, such 
as DTI and HARDI, should be enhanced such that fiber junctions are maintained, while reducing high frequency noise 
in the joint domain x S'^ of positions and orientations. Therefore we have considered diffusions and Hamilton- 
Jacobi equations on DW-MRI HARDI and DTI induced by fundamental stochastic processes on x S'^ embedded 
in the group manifold SE{3) of 31) -rigid body motions. In order to achieve rotation and translation co variance, 
processing must be left-invariant. We have classified all left-invariant linear and morphological scale spaces on DW- 
MRI. This classifications yield novel approaches, that in contrast to the existing methods, employ fiber extensions 
simultaneously over positions and orientations, while preserving crossings and/or bifurcations. Application on DTI 
creates fiber-crossings by extrapolation, comparable to HARDI, allowing to reduce the number of scanning directions, 
whereas applying the same method to HARDI removes spurious non-aligned crossings. In order to sharpen the DW- 
MRI images, we apply morphological scale spaces (Hamilton- Jacobi equations) related to probabilistic cost-processes 
on R^ X 6*^ . Moreover, we provide pseudo-linear scale spaces combining morphological scale spaces and linear scale 
spaces in a single evolution. All evolutions can be implemented by stable, left-invariant finite difference methods, 
or by convolution with analytic kernels. For example, we have shown by extending existing results on R^ (and 
recently on the Heisenberg group to SE{3), cf.|53 1), that morphological R^ x 5^ -convolutions are the unique viscosity 
solutions of Hamilton- Jacobi equations on R^ y\ S'^ and moreover we have provided analytical approximations of the 
morphological kernels (Green's functions). Finally, we have extended anisotropic linear diffusions systems on R^ y\ S'^ 
to data-adaptive diffusion systems of Perona Malik type that are useful to avoid undesired interaction between more 
isotropic regions (at the neural fibers) and less isotropic regions (at the ventricles). 
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In future work we aim to extend our differential geometrical approach on x S'^ to the actual tracking of fibers. 
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A Metric and Lagrangians on (M^ x S'^, dA^, dA^) 

Recall from Figure [6] that our left-invariant diffusions on x S'^ take place on the contact/sub-Riemannian manifold 
{SE{3),dA^,dA'^, dA^). On this contact manifold {SE{3), dA^.dA^, dA^) we set the Lagrangian 

where r] > ^ and where D^^^D^^ > and where 7*(s) = (d^*|^^^^ jli^)) ^re the components w.r.t. the moving 
frame of reference attached to the curve, cf. Figure [6] For 7^ ^ 00 we have a homogeneous Lagrangian and a 
left-invariant semi-metric on (6'£^(3), d^^, d^^, d^^) is set by 

^(^1,^2) := d{g2^gi,e = (0,/)) = 



inf inf /' ^\dA^{j{s))\^ ds ..... 

hiMe{0}xSO{2)cSE{3) 7 = (x(-),J?n(-)) e c/°°((o,L),5s(3)), ^ ]! i^^jl^^^^y ^^^^^ 

7(0) = eh,2 , 7(1) = 92 91^1' ' ' ' 

I ^ I 'Y I 

By left-invariance of the Lagrangian and t he fact that = we may as well set 7(0) = e and 7(L) = 
92^ 91^1 = {92h2)~^9ihi in Eq. (115) so that ^{§1^02) = ^(^1^1, ^'2^2) iov all hi,h2 G {0} x S0{2) and 
thereby inducing a well-defined metric 

d''"((y,n),(y',n')) = c)((y,i?„),(y',i?„')) 
on X S'^. Along horizontal curves s ^ 7(5) = {x{s), R{s)) in SE{3) we have 11(5) = R{s)ez = x{s) and 

\\K{s)r = \{dA\^^,j{s))? + \{dA%^^,iis))\^ 



{dA%^^^,i{s)) = \\x{s)\\ = l, 



and thereby the metric reduces to 



dh-((y,n),(/,nO)= inf / a/ + ^ d. (116) 



x(-) G C°°((0,L),R"^), 
x(0) =0, i(0) =e,, 
x(L) = i?„^(y'-y), 



where 5, L > 0, and n{s) are respectively spatial arclength, total length, and curvature of the spatial part of the curve, 
o te tha i 

malizec 



No te tha t the distance ( |116| ) is inde ed in variant under the choice of Rn such that Eq. ([9]) holds, due to left-invariance 
of (|ll5|) and rotation co variance of ( 116 ) and {RnRe^,a)^ = Rl^,a^n- Furthermore, we note that torsion is implicitly 
penalized by simultaneous penalization of length and curvature. 
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B Viscosity Solutions of Hamilton- Jacobi Equations on (M^ x S^, dA^) 



After introducing a positive semi-definit^ left-invariant metric tensor G : x 5^ x T(R^ x 5^) x T(R^ x S*^) ^ 
given by Eq. (36), we considered the solutions (51 ) and (53 ) of the Hamilton- Jacobi equations on x 6'^: 



dW 
dt 



(y,n,t) = ±(i/(dVr(-,-,t)))(y,n) := i(G(,,„))-i (dVr(-, •,t)|y„ , dW^(-, ^^^^^ 
[ Vr(y,n,0) = f/(y,n) 

where H stands for the Hamiltonian. It is well-known f20l that Hamilton- Jacobi equations in general do not have a 
unique solution, unless extra requirements such as the practically reasonable viscosity requirement are added|^ 



Our claim was that the solutions ( [5T] ) and ( [53] ) are the unique 1201 viscosity solutions of the Hamilton- Jacobi equations 
( [T36| ) and ([52]). 



Definition 8 Suppose that the Hamiltonian is convex and H{p) oo as p ^ oo. Then a viscosi ty so lution is a 
bounded and continuous (not necessarily differentiable) weak solution W : (R^ x 5^) x R+ M of (117) such that 



1. for any smooth function V 
one has 



2. for any smooth function V 
one has 



X S*^) X R+ R such that W — V attains a local maximum at (Jq^ ^o, ^o) 
dV, 



dt 



(yo.no, to) T {H{dV{; t)))(jo,wo) < . 



(118) 



y\ S'^) X R+ R such that W — V attains a local minimum at (Jq^ ^o, ^o) 
dV , 



dt 



(yo.no, to) T {H{dV{; t)))(jo,wo) > . 



(119) 



We will first provide a quick review of Hamilton- Jacobi theory, the Hopf-Lax formula and erosions on respectively 
R^ and on the Heisenberg group (5'E'(3))o that we obtain by contraction from SE{3) and which serves as a local 
approximation of SE{3), as explained in Section 10.1 



Definition 9 Let X be a normed space, then the Legendre-Fenchel transform L ^ on X is given by 

idxL){x) = sup {(x, y) - L{y)} 
yex 

where L : X ^ R U oo and (x, y) = x{y) and x G X*. In particular if X = one gets 

{dRr.L){x) = sup {x-y - L(y)} 

and in case of the Lie-algebra C{SE{3)) of left-invariant vector fields on SE{3) we get 

6 6 6 6 

{dciSEi3))L){^PidA') = sup {{^pidA\^c^Aj) - Li^c'Ai)} = {{dRd){pi, . . . ,pe) 

i=l j2 c-Ai 

i = l 



withl{c\...,c^) = cUi). 

^It need not be strictly positive definite, but the Hormander condition f48l should not be violated. 

^Consider for example the initial value problem ^ {x, t) = ( ^ (x, t))"^ and [/(x, 0) = 0, x G M, t > 0, with continuous (non-sensible) 

solution U (x, t) = \x\ — t if t > \x\ and zero else. This "nonsense" solution is not a viscosity solution, namely consider (f){x, t) = e~*e~^^ and 
= — 1 in the Definition of viscosity solutions according to Crandall l20l . 
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Lemma 1 Let t] > \. The Legendre-Fenchel transform of the function 

3 {c\c^,c\c^) ^ 



27^ V 
is given by 

9 (Pi,P2,P4,P5) ^ ^ {D'\{Pif + ip'?) + D^^ip^f + {p^fY e R 

Proof It is well-known that the Fenchel transform of c ^ ^11^11" equals p ^ ^||p||^ with \ ^ \ = 1, furthermore 
^oVx = Vx-i o 5^ where Vx denotes the scaling operator given hyVxf{x) = /(A~^x) from which the result follows 
by setting b = 2r], a = i^pi and respectively A = \fD^ ^ \fD^. □ 

Corollary 1 Let r] > ^, t > 0, then on the Hopf-Lax formula is given by 

w{x,t) := inf inf { / Cr,{^{s)) ds ^ u(y)} = minHjCr, (li^-y)] (120) 

7(t) =x 

7(0) =y 

with w{-^0) = u a given Lipschitz continuous function on R^, with 



(121) 



and w{x^ t) given by {120) is the unique viscosity solution of the Hamilton- Jacobi-Bellman system 
[ w{x^ 0) — u{x) 



with X = {x\x\x\ X ) G R , t > 0. The morphological Green' s function of this PDE is given by 

In the limiting case 77 ^ oo the Lagrangian is homo gene ou^^and 

w{x,t) = unn{u(y) + {{x' - y'f + {x^ - y^) + ^ {{x^ - y' f + {x^ - y^)^)^ } 

and we arrive at the time -independent Hamilton- J acobi equation 

In the limiting case r] ^ ^ the Hamiltonian is homogeneous and the Hopf-Lax formula reads 

w{x^t) = mm{k^ ' \x—y)^u\y)] 



where the fiat morphological erosion kernel is given by 



D",i,«,i/2,+(^) ^ I if^^^^^^ + ^^^^ < 

^ 1 00 else . 



The corresponding Hamiltonian 1-L{p) = is homogeneous as well, but now with respect to := | rch:3l l67l . 
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Proof Directly follows by Lemma [T] and (351 ch:3(Thm 4), ch:10 (Thm 1, Thm 3)]. For the special case of a 
homogeneous Lagrangian (r] oo), see \6T, Ch:3]. 

Remark 5 We parameterized elements in by means ofy = (y^, ^y^^y^) since we want to stress the analogy to 
Hamilton- J ac obi equations on y\ S'^. 



Lemma 2 (semigroup property of erosion operator on R^ y\ S'^) 

Let T] > \ and let D^^ > 0, D^^ > 0, t > 0. Let the morphological Green' s function be defined a^^ 



inf 

7 = (xi-), /?(•)) e c°°((o 

7(0) = (0,1 = Re^), 7(*) = (y, Rn), 

{ dA^\ , 7) = ( dA^\ , 7) = 
I7 I7 

inf / £^(7(5), 7(5)) d5 

7 = (^(O, -R(-)) e c°°((o, ^-i(t)), S£;(3)), q 

7(0) = (0,1 = Re^), 7(^~^(t)) = (j, Rn), 
{ dA^\ , 7) = ( , 7) = 



(122) 



with 



lip) 



4(7(s), 7(s)) = £,(( |^(^^ , j{s)), ...,{ dA\^^ , j{s))) 



217-1 

2v ■ 



(123) 



(124) 



where we recall Eq. {121 ) and with R^ x S'^- ''erosion arclength" given by 

P{r) = J VG-rir) (7(r), i{f )) df = J l^l(d^^l,(.) , 7(r))P df 

Y iG{l,2,4,5} 

with D^^ — D'^'^ and D^^ = D^^ and with left-invariant metric tensor G(-, •) = lim (>C„(-, •))^ and where iIj{s) = p 

77^00 

given by 

(125) 



V'(s)= / ^/k^~s) + /3^ ds. 



Then we have the following identity 



QU) = k^ 

■2 



eu 



(126) 



for all r G [0, t] and all Lipschitz continuous functions U : R^ x 5^ ^ R, where we recall that our erosion operator 
is defined in Eq. (39). 



Remark 6 We apply the following convention regarding curves within SE{3): Eor concise notation we shall write 
7(p) instead of J {s{p)). When writing j{p) we mean ^l(}l^~^(jp))y "^hen writing 7(5) we mean ^7(s) / 



Proof First of all we note that D^^ = D'^'^ and D"^^ = D^^ are necessary requirements for a well-defined morpholog- 
ical kernel '^'^ : R^ x 5*^ ^ R+ and in particular (as 7^ ^ 00) a well-defined metric on the group quotient 
R3 X 52 = (R3 X 5O(3))/({0} X 50(2)). 



Secondly, we will explain and derive the identities Eq. 
the second entry so that 



122) and Eq.( 123 ). The limit lim £„(•,•) is homogeneous in 
' II 77^00 



1 and lim k 

77^00 



(y,n) = t, (127) 



^^We assume that j{t) is chosen such that a minimizer exists, which is the case if j(t) is close enough to the unity (0, /). However, one can 
show that for 77 — ^ co and 7(f) far away from the unity a minimizer may not exist, similar to the 2D-case {12 \. 
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so t > is the total length of the curve in the sub-Riemannian manifold (R^ x S*^, d^^) and where p is the arc-length 
parameter in (R^ x S*^, d^^) since connecti ng cu rves over which is optimized are not allowed to use the ^3 -direction 

in the tangent bundle. Furthermore, by Eq. ( 127) we have that ^£oo(7(p), "^"^ = ^2^' which explains 

Eq.( 123 ). Regarding Eq.( 122 ) it is relevant to keep t > fixed during the curve optimization in ( 122 ) for ^ < r] < 00, 
otherwise the minimizing curve may not exist, and direct computation yields 

/£,(7(P),7(P)) (af)"^ dp = /£,(7(p), 

^ ^ 

where we used the homogeneity of Crj : 



dr^dp\ ds — 
dp ds ) dp^^ ~ 
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dp ^ ^ dp ds 



Thirdly, we recall that functions U : R^ x 5^ ^ R are related to functions U : SE{3) R with invariance property 

0^geSEi3)yh^io,R,^,^)eSEi3)U{gh) = /7(^)) by means of 



with g = (x, R) e SE{3),y G R^ 



U{g) = U{x,R) = U{x,Re,), 
= (0, 0, 1)^. We apply this identification by setting 
W{g,t) := Jg^^^{h{q-'g) + U{q)} 



for all g G SE{3) and t > 0. Consequently, we may rewrite our Eq. ( [126| ) as 

W{g,t) 



(128) 



where we use short notation /ct(y, R) := 



Fourthly, we apply basically the same approach as in [53, Thm.l] where we replace the Heisenberg group by the 
3D-Euclidean motion group SE{3) where in this lemma we carefully omit scaling properties of the Lagrangian that 
do hold for the Heisenberg group, but not on the SE{3) (even though one has the contraction \imq^o{SE{3))q = 
{SE {3))o which serves as a local nilpotent approximation, |68|). We shall return to these local appro ximations later 



when we derive approximations of 
In general one has 



(y, n). For, now let us proceed with proving Eq. ( 128 ). 



kt{v ^g)<kr{v ^q)^kt-r{q ^g) 



(129) 



for all g^h^v G SE{3) and all r G [0, t] as the concatenation of two optimal curves yields an admissible curve over 
which is optimized in the left hand side. Here equality is obtained if is on the minimizing curve between e and 
v~^g, i.e. by left-invariance g is on the minimizing cur ve be twee n e an d v~^g. Now due to continuity and convexity 
of (c^, c^, c^, c^) ^ Crj{c^^c^^ c^, c^) the infimum in ( 128 ) and (122) is actually a minimum and therefore we can 
choose V G SE{3) such that 

W{q,T) = kr{v-^q)^U{v). 



Then by Eq. ( 129 ) one has 

W{g,t) <kt{v-'g)^U{v) 

<~kr{v-^q) ^k-r{q-^g) ^U{v) 
= W{q,r)^k-r{q-'g) 
for all q G SE{3) and thereby (by taking the infimum over all q) we obtain 



QU)>kr^''^^^'^^eu. 



(130) 
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(fci>",i>",.,+ © C7) < kr'"' e C/. Let 



So in order to prove ( 126 ) it remains to be shown that k 
w e SE{3) such that 

W{g,t)= min kt{v~^g) ^ U{v) = hiuj-^g) ^ U{w) , 

veSE{3) 

now take g as a point along the minimizing curve between g and w where the vertical y\ S'^ arc-length equals p = r. 
Then we have 

W{g, t) = kr{w-^q) + k-r{q-'g) + U{w) > k-Aq-'g) + W{q, r) ^ 

> min \kt-T(r~^Q) -\- W(r.r)} 

reSE{3) 



from which the result follows. 



□ 



Lemma 3 Let {Ai}^^-^ be a basis of the Lie-algebra of tangent vectors at the unity element and let [Aj] be the 
corresponding left -invariant vector fields and let {d^*} be the corresponding dual vector fields, e.g. Eq. {33). Then 
the exponential curves in SE{3) 

6 

j{h) = {x{h), R{h)) = exp{h ^ c'Ai) (131) 

are auto-parallel w.r.t. Cartan connection. Their tangent vectors have constant components w.r.t. the attached left- 
invariant moving frame of reference, i.e. 



for i = 1, . . . ,6 and all /i G R. The spatial arclength of the spatial part x of the curve is given by 

s{h) = h. 



(132) 



.i\2 



Let = = 0, then the sub-Riemannian {SE{3)^ dA^^ dA^) -arclength ( if = c^, and = this corresponds 



to R X S - ''erosion arclength", Eq.{55)) is given by 



p(h) = h 



\ 



^G{1,2,4,5} 



,i|2 



Letc^ = = = 0, then the sub-Riemannian {SE{3), dA^^ dA^^ dA^) -arclength is given by q(h) = h I ^ 



iG{3,4,5} 



Proof By definition, 7 is given by Eq. (131) the unique 1 -parameter subgroup of SE{3) whose tangent vector at 



the unity element e equals ^ cM^. It is thereby the integral curve of the corresponding left-invariant vector field 

6 6 6 

given by (Lg)^ c'^Ai = ^ c'^{Lg)^Ai = c'^ Ai\g, obtained by push- forward of the left-multiplication, yielding 

i=l i=l i=l 



= from which Eq.(132 



i=l 



follows. In the Cartan connection on the tangent bundle T{SE{3)), the 



antisymmetric structure constants serve as Christoffel symbols, cf . 1311 and we have 

V^7 = 0^V,=i,...,6 : 7' + ^cife7^7'==f =0^V,=i,...,6 : f := {dA\ , j) = c\ 

For spatial arclength we have ||x(s)|| = 1 <^ ||x^(/i)||^ = (c^)^- Regarding arclength in the sub-Riemannian 

^e{l,2,3} 

manifold {SE{3), dA^^ dA^) we have by Eq.( 132) that along a horizontal exponential curve 
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y 2e{i,2,4,5} 

gously. 



and the other sub-Riemannian case can be derived analo- 

e{l,2,4,5} 

□ 



Let ' ' respectively denote the viscosity solution of Eq. (136), Eq. (52), with initial condition ±5c, rj > 



\,D^^,D'^^ > 0. This kernel describes the propagation of "balls" in x SE{3)/{{0} x S0{2)) centered 

around unity element e = (0, e^) = {(0, Re^,a) | G [0, 27r)} G ><\ S'^. Next we derive a few estimates for this 
morphological kernel in order to 



1. motivate our analytic approximations in our erosion algorithms described in Section [7] 

2. serve as ingredients in Theorem]?] which is the main result of this section. 

Lemma 4 Let r] > ^, D^^^ D^^ > 0, h > 0. The morphological kernel '^'^ satisfies the following estimate 



iG{l,2,4,5} 



expj h Yl 

e{l,2,4,5} 



(133) 



iG{l,2,4,5} 



with Ai = Ai\^ e Te{SE{3)), c' e R, C > 1 and \ ^'M = / E w^W? ^ith D^^ = D'^ 

iG{l,2,4,5} yiG{l,2,4,5} 



and D^^ = D^^ and where 



exp \ h Y 

ze{l,2,4,5} 



denotes the left coset in W x = SE{3)/{{0} x S0{2)) 



associated to the group element exp \ h Y ) e SE{3). 

ze{l,2,4,5} 



Proof We first derive an explicit upper-bound. As the smooth exponential curve 



[0, h]3h^ 



exp ih ^2 ^^^i 

iG{l,2,4,5} 



G X 



is horizontal and connects [el and 



exp \h Y 

ze{l,2,4,5} 



7-)ll r)33 1 

we have by the definition of kj; ' '^'^ in Eq.(|122|) that 



exp h Yj ^^^i 

\ iG{l,2,4,5} 



< — C72r?-i I 2^ c'Ai\'^^- 



h 



2ri 



iG{l,2,4,5} 



with 



G{1,2,4,5} 



where we applied Lemma 



and used ^ = ^jj^ along the smooth horizontal exponential 



curve. The lower bound follows by the theory of weighted subcoercive operators on Lie groups, (681 ch: 1,6] where one 
should consider the special case U = 1Z, G = SE{3) and algebraic basis {^i, ^2, ^4, A^} as erosion/dilation only 
takes place in the directions Ai, A2, A4 and A5 (i.e. G-orthogonal to ^3). As a result we must assign the following 
weights to the Lie-algebra elements wi = W2 = = W4 = 1 and ws = wq = 2. □ 
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Clearly, not every element (y, R^^ ^ SE{3) is reached by an exponential curve of the type h e 



{1,2,4,5} ' 



In fact, by the Campbell-Baker-Hausdorff formula and the commutator table (29 ) one has 



= exp{h{y^Ai + ^^2)} exp{/i7A4} expihpA^} 

= exp{h{y^Ai + y^A2)}exp{hjA^ + f3hA^ + ^(3h^[A^,A^] + 0{h^)} 

= eMHy'AiWA2) + h^A^ + (3hA, + ^h\y' p-y^^)As + ^h\(3^)Ae + 0(/i3)} 

6 

= exp{ c'Aih^^ + 

i=l 

(134) 

and consequently, we can extend the estimate in Lemma |4j 

Lemma 5 Let be a compact set around the unity element e = {O^Cz) = [{0^1)] = {{O^Re^^a) \ <^ G [0,27r)} 
within x S'^, then there exists an e > and C > 1 such that for all h < e one has the following uniform estimate 
on compact sets 



for all {y^n{f3^'y) = [expj^ c^Aih'^^}] G ft^, withwi = W2 = = W4 = landws = wq = 2 and := c*(j,a 



0, /3, 7) = c*(j, a = 0, /3, 7), given by Eq. { 75) and { 72). 



Proof Regarding the assignment of the weights to the Lie-algebra elements we refer to Eq. (134). Regarding the 
logarithmic mapping in SE(^) we recall Section 9.1 Akin to the metric (162), the formula (122) is well-defined 
on the partition -a S'^ := SE{3)/{{0} x S0{2)) of left cosets. However as we have the restriction (d^^, 7) = 
(da, 7) = (da, 7) = we must apply a consistent cross-section in x S'^. For convenience, we take the unique 
element from the left cosets with a = in the Euler-angle parametrization in both the initial point (0, R^^ = /) and 
endpoint (y, R^ = R^^^^R^^ ^) of the curve. For the rest the proof is the same as the proof of Lemmajij □. 



Corollary 2 Let r] > ^, D^^ > 0^ D^^ > 0. wi = W2 = w/^ = wr^ = 1, w^, = wq = 2. For the morphological 
erosion (-\-) and dilation kernel (-) onM.^ y\ S'^ one can use the following asymptotical formula 



\i=l 



2r] 



(135) 



for sufficiently small time t > 0. 



Lemma 6 Let G R, i = 1, 2, 4, 5 and t] > \ then 



lim 



exp \h 

V iG{l,2,4,5} . 



27^-1 



h 



E 



277 ' ^-^ \/D^^ 
' 2,4,5} 



d , , 277 277 —1 

--AA^r,-! - 



27? 



/ ^ ]Jii 

zG{l,2,4,5} 



,i|2 



Proof Consider Eq. (122), then to each C -curve with (d^^|^ ,7) = (d^^|^ ,7) =0, we have the following 
-functions 



[0,/7] 3p^f{p) := (d^*|^^^^,7(p)) G 
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for i G {1, 2, 4, 5}. By applying a first order Taylor- approximation around 5 = we obtain 

= lim/i-i / Cr,{^j\0) + O{s),^^H0) + O{s),^^\0) + O{s),^jH0) + O{s))ds 
and T] > ^ implies that this is of the order 

' \iG{l,2,4,5} 

where ?/^~^ (0) = and along a minimizing curve we have finite curvature (as we will see in Appendix [g]) so ?/^~^ (/i) = 
0{h) and the first term vanishes as | and the result follows by definition of the morphological kernel, Eq.( |122| ). □ 

Theorem 4 The viscosity solutions of 

-(y^n^t) T I (g^,^,) (diy(-, .t)!^,, , dW{-, J)' = 
W(y,n,{)) = U(y,n) 



dW , 
dt 



(136) 



with G(j^n) = 9ii ^^^\y n ^ ^^^\y „ ^^^^ ~ ^^^^ ~ ~ ^ ^ respectively given by (-\- 

case) left -invariant erosion 



(fcr ■''■"eK3,5.C/)(j,«)= inf U(y',n') + kr -/), i?» 



(137) 



and (— case) left-invariant dilation 

(fcf (Bu.^s^ U)(y,n) = sup [fcf "■^"•''•-(i?„^(j + U(y',n')\ (138) 

Proof The proof consists of two parts, first we must show that they are indeed solutions and then we show that they 
are viscosity solutions. We will only consider the erosion case since the dilation case can be treated analogously. 

part I. This part consists of two subparts. In part la we will show that if we set 

W{y,n,t) := (fcf"'^"''''+ Qm^s^ U){y,n) 
with morphological kernel given by ( |122| ) then 

^{y, n,t) T 2^ (G^y^n) {dW{; , dW{; ; t)Qy < . (139) 

Subsequently, in part lb we show that 

^(y, n,t) T 2^ (G^y^n) {dW{; ■,t)|^_„ , dW{; ; t)Qy > . (140) 

so that part I is finished. 



part la Again we resort to evolutions on the full group SE{3) by setting VK(y, t) = VK(y, R^z^t) for all y G 
R e S0{3), t > 0. Then for all A = c'Ai one has by LemmaSthat 

iG{l,2,4,5} y 

Wige^'^.t^h) = min ^^'''^'''"^'+(^-1^/6^^) + iy(v, t) 

v^SE{3) 

<^f' •^"•''•+(e''^) + H^(5,t) 
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for all g = {y^R) and consequently, one has 

h - h 

and thereby, by Lemma [6] and the construction of the left-invariant vector fields by means of the derivative of 
the right regular representation Ai = d1Z{Ai) one gets 

^^jf ^^ +c' A\^W{g,t)<C^{c) (141) 

for all c = {c^ , , , c^) G M^. So by subtracting the Lagrangian and taking the maximum over all c we apply 
the Fenchel transform of the Lagrangian which yields the Hamiltonian, i.e. 



I iG{l,2,4,5} 
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dWig,t) 1 J ^ \Ai\gWig,t)\ 
dt ^ 277 



^{1,2,4,5} J 



part lb Let g* = (y*, i^n* ) ^ SE{3) be the minimizer in the erosion operator, i.e. 

W{g, t) = W{g\Q) + fcf 

Then we have 

-(w^(5,0)+fcf"'^''>''>+((5*)-i5e-^^)e+*^*)) (142) 

6 6 

where ^ = e = e^=i and g* = e = e^=^ G SE{3). If we now let ^ then we obtain 

= limi sup C,'^ + ,,,,, 
HO yvG5£;(3) (143) 

> lim ^ 

where we applied Lemma [2] and where we note that by the Campbell-Baker-Hausdorff formula 

g-f Ag+f A* ^ g-t(A-A*) + i^[A,A*]+0(/i3) 

SO that 

V^(^e-^(^-^*\t - h) = I^(^e-^^e+^^*,t - /i) + 0{h^) 
for almost every ^ G x S*^ (where W{-^t — h) is differentiable). 

Recall that the left-invariant vector fields are obtained by the derivative of the right-regular representation Ai = 
d7l{Ai), Ai = Ai\ J) and as a result the limit in the left-hand side of ( 143 ) equals 

iG{l,2,4,5} ^ ^ 
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Furthermore, we can assume without loss of generahty that = = (c*)^ = (c*)^ = 0, i.e. 

A - A* G span{Ai, ^2, ^4,^5}. (144) 
To this end we recall both x 6'^ = 6'^(3)/({0} x S0{2)) and W{g, t) = W{y, Rez.t), g e SE{3) so that 



for all a, /3 G M and by the Campbell-Baker-Hausdorff formula one can always find a, /3 G M such that 

log (e^e"^^) - log (e^*e^^^) C span{Ai, A2, A4, A5}. 



Assumption ( |144| ) allows us to compute the limit in the final right-hand side of inequality ( |143| ) by means of 
Lemma |6] and we obtain 



2,4,5} ^ ^ 

(A-A*)x 



> '-^(-^r^') = (^^^, ^^^^^), ^^^^) . 



from which we conclude : 



iG{l,2,4,5} 
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2,4,5} V / ^ V 



\ iG{l,2,4,5} 



c — (c j c — (c ) \ c — (c 



from which the result ( |14Q| ) follows. 



/?ar^ //. Next we verify that erosion ( 138 ) with the Green's function k 



(yo, no, to) e (M^ y\S'^)x where W — V attains a local maximum. The other case ( 1 19 ) can be shown analogously. 



Similarly one can show that dilation (39) with the kernel /c^ ' ' 
Eq. ([52]). 



indeed satisfies ( 1 1 8 ) at a location 



is the viscosity solution of 



First of all it follows by left-invariance (and the fact that SE{3) acts transitively on x S'^ by means of Eq. ( 18 )) that 
without loss of generality we can restrict ourselves to the case (yg, no) = (0, e^) and furthermore, by the semigroup 
property Lemma[2j we can, again without loss of generality, restrict ourselves to the case to = 0. Since W — V attains 
a maximum in (0, e^, 0), there exists a small open set Q around (0, e^, 0) in (R^ x 6*^) x R where 



W{y, n, t) - F(y, n, t) < W{y, n, 0) - F(y, n, 0) 
Furthermore we have by Eq. ( |137| ) that 

W{y,n,t) > W^(0,e,,0) + fcf"'^"''''+(y,n) , 



for all (y, n) G xi S^. Combining the estimates ( 145) with ( 146 ( yields 



V{0, e,, 0) + V{y, n, t) < W{y, n, t) - W{0, e„ 0) < k'^ 



(y,n) 



(145) 



(146) 



(147) 



locally around (0, e^, 0), that is within f2. Or equivalently for the corresponding function on SE{3) given by V{x, R) = 
V{x, Rez ) we have 



-y(0,J,0) + y(y,i?„,t) ^ fcf •''•+(y,n) 



< 



(148) 
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for t e [0,ei),d5£;(3)(y,i?n, (0,/)) < 62, for some ei > 0, 62 > 0. Now set h = {y,Rn) = e^^ with A = 
^ d'Ai and take the limit | then we have 

iG{l,2,4,5} 

^(0,/,0)+ ^ cM7^(Ai)F(0,/,0)<^^ I ^ cMil^ 

iG{l,2,4,5} ^ iG{l,2,4,5} 

where we appHed Lemma [6] (which tells us that the constant C in Lemma|4]can be set to 1 in the limiting case h | 0). 
Now the left-invariant vector fields are given by the derivative of the right regular representation Ai = d1Z{Ai) with 
Ai = Ai so we see that 

^(0,/,0) + sup {cMil/(0,/,0)-^| E cMi|5|^}<0 

.eu^.)^'-"' '^^'''''''^ (149) 

f (0,/,0) +{^^^^sEi3)){A^ '-^\A\^^)){dV{;0)) <0 

5 

with A= ^'^i ^ Te{SE{3) and dV'(-, t) = AiV{-, t)dA' e {Te{SE{3))y. Now by the final remark 

iG{l,2,4,5} i=l 

in Definition |§] and Lenmia[T]we see that 

f (0,7,0) + (^?R4((c^c^c^c5) ^ 

(Eie{i,2,4,5} Icf ) 2;^)) {AiV{fi, 1, 0), (^21/(0, /, 0), {A^ViO, 1, 0), (^5^(0, /, 0)) < (150) 

^ §(0, /, 0) + 2^ ((AiF(0, /, 0))2 + (^2^(0, /, 0))2 + (^4^(0, /, 0))2 + (^5^(0, /, O))^)'' < 

As a result we conclude that the erosions Eq. ( |137| ) with morphological kernel ( |122| ) are the (unique) viscosity solutions 
of the Hamilton- Jacobi equations ( 136 ) on the contact manifold (R^ x S*^, d^'^). □ 



C The Hamilton- Jacobi equation on (M^ x 5"^, d^^) and the propagation of 
geodesically equidistant surfaces in x 

Recall from Figure [6]that our left-invariant erosions on x take place on the contact/sub-Riemannian manifold 

iSE{S), dA^, d^6). On iSE{S), dA^, dA'^) we set the Lagrangian (for r/ > i) 

c,h{p)Mp)) = ^ [ + ^44 ) ■[ds) ' ^151) 

expressed in the (S'^(3), d^^, d^^) arc-length parameter p given by 



Pir) =!J E i^|(d^'|^(,),7(T))Pdf. (152) 

y iG{l,2,4,5} 



yielding (when rj oo) the well-defined metric on (R^ x /S^, d^^) given by 

d((y,n), (y', n')) = d(((i?„')^(y " (i?n')^n), (0, e.)) , 



Pmax 

^^((y,n),(0,e^)) = inf J Coo{l{p)A{p))^P' (153) 

7(0) = (0, ), 7(Pmax) = (y, J^n), 
< d^3| , ^) = ( dA^\ , 7) = 

I7 I7 

Let us return to the general Lagrangian minimization problem on (R^ x S'^^dA^) for 7^ > ^ 

Pmax — ^ / \ 

inf f >C„(7(p),^(p)) • (^) dp. 
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for some given go G SE{3) and where the total arc-length p^ax in x 5*^ is fixedj^ say pmax = t > 0. Along 
an optimizing curve 7* = (x*, i?*) ("geodesic") in SE{3) with corresponding curve s ^ 7*(5) = (x*(5),n*(5) := 
R*{s)e,) G X we have 



|E,(7*(t),t) = /:,(7*(t),7*W) • 



where we recall that iIj{s) = p and s = i/j ^ (p), Eq. ( 125 ). Now consider a family of surfaces associated to a smooth 
function W : x 5'^ x R+ ^ R+ given by 

St := {(y,n) G X 52 I I^(y,n,t) = Wo} 

parameterized by t > 0, where Wo > is some positive constant. Such a family of surfaces is called geodesically 
equidistant, 1. 67 J , if 



dw , 

dt 



(dW^(-,*)lr(t),rW) = (Ve£,(7*W,c)-c)|^^..(^^. 

The chain-law now gives 



(155) 



dt 



where according to the 2nd equality in ( |155| ) we may rewrite the righthand side as 

sup ^ {-£,(7*(t),c) + (dVr(.,i)|^.(,),c)} 

ie{l,2,4,5} 

which equals the Legendre-Fenchel transform of >Cr7(7*(^), •) on T^*^i^{SE{3)) that can be computed by the results 
in Appendix |b] recall Def. [9] and Lemma[l]so that we obtain the Hamilton- Jacobi equation on R^ x S*^: 



dW I 
dt 



(7*W,^) = ^ (G-\,)(diy(7*(t),t),diy(7*(t),t)))' 
-^{Y{t).t) = ^^[D'HA,\^.^,^W{r{t).t))'^D''{^^^^ 
+I)^^(A|^.(,)I^(7*W,^))'+^'nA|,*(,)m7*W,^))^" ^ 



along the characteristic curves. So we conclude that iso-contours of the solutions of W in Eq. ( |52|) and Eq. ( |136| ) are 
(at least) locally geodesically equidistant. This allows us to use the Hamilton- Jacobi equations ( |l36| ) for wavefront 
propagation methods |[62| |56l for finding geodesies in R^ x S*^ . We will consider this in future work. 



D Asymptotical Expansions around the Origin of the /c-step Time-integrated 
Heat Kernel on x 

Recall from Section [T0]2] the Gaussian estimates for the heat-kernels (contour enhancement); 

iff ''^ (y,n) = 



167r2(D33)2(£,44)2i4 



Otherwise for fixed 77 there does not exist a minimizer. 
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with 1^1 given by Eq. ( 87 ). Then we obtain the following relation (with again short notation g = (y, R^) G SE{3)): 



R 



■X,k 



(y,n) := {{-D^^A^f - D^^A^f - D^\A,f + XI)-'^5,) (y,n) 



(y,n) r{t; k,\) dt 







dt 





'^\g\'^-^K,_k{\9\y/X) 



— ^2(^33)2(^44)2 (^_i) 

with A52 = A\ and now we have the following asymptotical formula for the Bessel functions 



- log(2:/2) - JEULER if = 



^ i(|H-l)!(f)-'- if^GZ,^^0. 
with Euler's constant 'Jeuler, which holds for < z << 1, so that for \g\ = |(y, Rn)\ « 1 we have 

|(;|-|fe-4|+fc-4A-^ 



R 



(y,n) 



7r2(D33)2(2:)44)2 | 
^2(^33)^(^44)2 (- log ■ 



\ 2~^A^ 
-2 lEULER) 3! 



if 7^ 4, 
if A: = 4 



Likewise in the contour completion kernel, recall ( 82 ), we get rid of the singularity at the origin in 

Rx,k = Rx,l *R3^52 ^A,l 

by setting k > dim(R^ x S'^) = 5. In other words to avoid a singularity in the Green's function one must use at least 
k = 5 iteration steps in the contour-enhancement process. 



E Conditions on a Left-invariant Metric Tensor on xi and the Choice 

of D^^ and g^^ 

In this section we derive sufficient conditions on a metric tensor ^ to be both left-invariant and well-defined on the 
quotient R3 x = SE{3)/{{0} x S0{2)). 

Definition 10 A metric tensor G : SE{3) x T{SE{3)) x T{SE{3)) C 

• is left-invariant iff 

Ggg{{Lg),Xg, {Lg),Y,) = G,(X„ Y,) (156) 

for all g,qe SE{3) and allX,Y e T{SE{3)). 

• provides a well-defined metric tensor on xi = SE{3)/{{0} x 50(2)) iff 

Ggh{{Rh)*Xg, {Rh).Yg) = Gg{Xg, Yg) (157) 

for all g G SE{3) and all h G ({0} x 50(2)) and all X, F G T{SE{3)), in which case the corresponding 
metric tensor on the quotient is then given by 

G(y,n){^^' A|(,,„) '^.■I(y,»)) ■■= G(y,R„){Ec' ■^Lr.)'^'^' M(,,rJ ^^^^^ 



62 



where vector fields are described by the differential operators on C^(M^ x S'^); 



U{y-\-hRnej ,n) — U(y—hRnej ,n) 



(^3+j|(^,„)^)(j,w) = lim ^ ^ ^ , J = 1,2,3, 



where Rej^h denotes the counter-clockwise rotation around axis ej by angle h, with ei — (1,0 ,0)^, e2 — 
(0, 1,0)^, ^3 = (0,0, 1)^. Note that the choice of Rn satisfying Eq. ([p]) does not affect Eq. {158) because of 

Eq.mn. 



Remark 7 In the other sections in this article, for the sake of simplicity, we do not distinguish between X and X. For 
example we write both As | U and As \(^y r^) U. 



Set A = {Ai^ ... ,^6)^ as a column vector of left-invariant vector fields on SE{3), then one has the following 
identity 



with Za G SO (6) is given by Eq. ( [46| ), which is straightforwardly verified by Eq. (31) where we computed the 
left-invariant vector fields explicitly in Euler- angles. As the functions U are invariant under right-multiplication with 
elements (0, R^J G ({0} x S0{2)) this yields the following identity 

(A U){gh)=t{Z^)){A, U){g), 



gh 



for all h = (0, Re,, a), where Z^ = Re,,a © Re,,a G S0{6), with i?e,,a ^ S0{3) is given by (|46]). Recall that by 
right multiplication with (0, Re,, a) one takes a different section in the partition of left-cosets in x 5^, boiling down 
to a rotation over a simultaneously in the two grey planes depicted in Figure [6] 



Theorem 5 A metric tensor G(^y^ri) {^^ x x T(R^ x S^) x T{R^ x S^) given by 



^iy,n) = Yl dijiy.Rn) dA' dA^ 
with QQj = QjQ = 0, j = 1, . . . , 6 well-defined and left-invariant iff 



which is satisfied iff 



^cx [9ij] = [9ij] , for all a G [0, 27r), and 
Qijiy^ Rn) = dij cire constant 

[gij] = diag{gn, gn, gss, 944, 944,0}. 



(159) 



(160) 



proof The left-invariant vector fields satisfy Ai = (Lg)^ Ai , since {Lg)^{Lq)^ = {Lgq)^. So the left-invariance 
requirement reduces to 

Lg*9ij = 9ij 

for all g G 5'^(3). Now 5'^(3) acts transitively onto SE{3) (and onto x 5^) so gij must be constant. Regarding 
( |157| ) we note that 



GghiiRh). E A , {Rh). E Aj 

i=l 9 j=l 

6 6 



6 

I 

i,j,k,l=l 



gh 



i,j,k,l=l 



i=l 



6 

E 
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for all (c\ . . . c^) e and all a e [0, 27r) so that 

VaG[0,27r) • ^a[9ij]^a = [9ij]' 

The final results follows by Schur's lemma on SO (2) (operators commuting with irreducible representations are 
multiples of the identity). □ 



The actual (horizontal) metric (induced by the metric tensor given by Eq. (10) on the contact manifold (R^^ x 

5^, dA^.dA^) is given by 

d(y, n , nO := d{R^. (y - i^?n , 0, e,) (161) 

with (horizontal) modulus given by 

/ , I E 9^J {^^%(p) A{P)) (d^^'U(p) dp , (162) 

^ = (x(.),R(.)) e c°°ao,t),SE(3)), ^ \ inaf^ 4.^,1 ^^^^ ^^^^ ^ 

7(0) = (0, 7 = Ke^),7(t) = (y, f^n), ^ y2,Jti^,4,D| 



Remark 8 Mz^^ measuring the distance (Eq. ( [7(57| ) and {162)) between the cosets [gi] = giH and [^2] = 92H, with 
H = {0} X S0{2) and = {xi^ Rg^^^-Re ^jSiRe^^ai)^ i = 1,2 one first determines the elements from both co-sets 
with vanishing first Euler angle, i.e (x^, Re^^^-Re^^f^.) and then compute the distance w.r.t. these elements in the full 
group where the connecting curves are not allowed to use the ''illegal" Aq- directions. Note that this procedure leads 
to a well-defined metric on (R^ x S*^, d^^, d^^j despite taking the section. The condition (d^^|^ , 7) = avoids 
possible short-cuts via the ''illegal" Ae-direction. 



E.l Data adaption and left invariance 

In the erosion algorithms one can include adaptivity by making I)^^ depend on the local laplace-Beltrami-operator, 
recall Subsection |7.l[ whereas in the diffusion algorithms one can include adaptivity by replacing D^^ in AsD^^As 
by a data adaptive conductivity 

D^HU)=e-'-^,K>0. 

In both cases this does not correspond to making a left-invariant adaptive inverse metric tensor, as this would contradict 
Theorem[5] However, an adaptive erosion generator (7^ = 1): 

5 5 
^g'^{U){y,n)-{A^0Aj){dU,dU)= ^ g'^ {U){y,n) • {AU){AjU) 

i,j = l «,j = l 

is left-invariant iff g'^^ is left-invariant, that is g'^^ commutes with left regular action of SE{3) onto L2(R^ x S'^). 
Similarly an adaptive diffusion generator 

5 

U^Yl MD''{U)Aj{U)), 

is left-invariant iff D^^ is left-invariant. 
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F Putting the left-invariant vector fields and the diffusion generator in ma- 
trix form in case of linear interpolation 



We will derive the N^No x N^No matrix form of the forward approximation ( 100) of the left-invariant vector fields 
and subsequently for the (hypo)-elliptic left-invariant diffusion generator. Here N stands for the number of spatial 
pixels (where we assume a cubic domain) and where No denotes the number of discrete orientations on S'^ . Note that 
the matrix-form of the backward and central differences can be derived analogously. 

Recall that the forward approximations of the left-invariant vector fields are given by 

p = 1, 2, 3 and where / G {1, . . . , A/'^} enumerates the discrete orientations Ui G S'^ and where y = {y^ ^ y'^ ^ y^) G 
{1, . . . , N}^, where h denotes spatial step-size and ha denotes denotes the angular stepsize. 



We approximate the angular vector fields by 

ha 



(^,+3t/)(y,nO ^ ^ f-/7(y,nO + ^ M,{;^-^+^ /7(y,n,), 



i'=i / 
where the upper-index f stands for "forward" and with 

J^f,ha,P+S ^ I ^^^^^^^ ^^^3^ 

[ else, 
where Ap^i is the unique spherical triangle in the spherical triangularization containing / := Rn^R^^ h^tz. 

We approximate the spatial vector fields by 

y^' ,2/^' ,y^' 

with /-indexed discrete spatial kernel given by 

171=1 

with (jp^i)^ the m-th component of the vector ^ := hR^^p, p = 1,2,3, and with linear interpolation kernel 
Va :Z ^ [0, 1] given by 

r l-|a| if 6 = 
Va{b) = \ H{ab)\a\ if6G{-l,l} 
[ else 

with heavy side function H{u) while assuming \a\ < 1. 

So if we store ((/7(y, n/))/^|i 7Vo})yG{i,.-,A^}^ one long column vector u G R^^^° we can represent the (forward) 
angular left-invariant vector fields by the matrix 

. / A^3 \ 1 

a J a 
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with (g) the Kronecker product and with M^^l^^ 

jn>N^NoXN^No 



the matrix with entries (163) and where In^No ^ 



denotes the identity matrix. 
The (forward) spatial left-invariant vector fields is now represented by the block matrix 



a=i,...,No 



with M/'^'^ e R^^^^^ the matrix with entries kf'^[y'^ - y^\y'^ - y'^\y^ - y^']. 

Similarly, one defines central (index by c) a nd ba ckward differences (index by b). The matrix-representation of the 
generator of the (hypo) -elliptic diffusion, Eg. (102), is now given by Jj^s + J52 G NoxN No ^j^j^ 



Js^ := D^^ (^AiAl + AiAl). 



(164) 



G Solving (the Pfaffian system) for geodesies on {SE{3) , d^^ , dA^ , d^*^) 



Here we follow the same approach as in IST] App.A] where we minimized on the contact manifold {SE{2)^ d^sE(2))- 
We extend the manifold SE{3) into Z = SE{3) x (R+)2 x R+ x {T{SE{3))y on which SE{3) acts as follows 

T]g^{x,y,z,R){{y^,y^,y^,B!),H.,cr, (A\ A^ A^, A^, A^ A^)) = (^(2/\2/\2/^,i?0,^,cr,CoAdp(A\ A^ A^, A^, A^ A^)) 

where CoAd is the coadjoint action of SE(3) onto the dual tangent space {T{SE(3))Y , where = x = —n^Ai + 
/^^^2 + 0^3 and ds = crdt = ||x'(t) \\dt represent curvature and arc-length s of the spatial R^-part of curves in SE{3), 
as by definition on this extended contact manifold (Z, 6^ ^ . . .6^) the following Pfaffian forms vanish 



0^ 



= d^^ - crdt = , 
= d^"^ - K^adt = , 
d^^ - K'^adt , 
= d^^ = 



(165) 



The first five Pfaffian forms have to vanish in order to ensure that 7 : [0, L] SE{3) yields a horizontal curve 
s {x{s)^n{s)) in R^ x S'^ where the spatial tangent corresponds to point on the sphere, i.e. 



±{s) =n{s). 



(166) 



The 6th Pfaffian form is imposed to avoid a spurious correspondence between horizontal curves 7 = (x, i?) in the 
sub-Riemannian manifold {SE{3), d^^, d^^, d^^) to curves (x, n) in R^ x S'^ satisfying (166) via 



{x{s),n{s)) = {x{s),R{s)e,). 

To this end we recall that AeU = for all smooth U : SE{3) R given by U{x, R) = U{x, Re^) with : R^ x 
S'^ R. and (d^^ |^ , 7) = excludes spurious curves in SE{3) by means of 

{x{s),R{s)) ^ {x{s),R{s)R^^^a^s)) 
with a 7^ 0, that would correspond to the same horizontal curves in R^ x 5^, see Figure 22 Our goal is to solve for 
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f (y, n) 



(0,e^) 



SE{3) 





A — F) 

71 










a 









(0,i?e 



7i ~ 72 , 



Figure 22: An (optimal) curve in x 5^ corresponds to a class of equivalent curves in SE(^). For convenience we 
impose the extra constraint ( | ^ ^ 7) = to avoid spurious torsion and curvature. In the figure we see that 71 and 
72 are equivalent, i.e. 71(5) = ^2{s)h{s) with h{s) = (0, Re^,a{s))^ but only 71 satisfies the convenient constraint. As 
the constraint = does not affect the final geodesic in x 5^, we can set the corresponding Lagrange-multiplier 
Ae = in Eq. ([173]). 



the minimizing curves of the optimization problem dllSk on {SE{3)^ d^^, d^^, d^^"^ 



^{91^92) '= ^{92 ^9i^e = (0,/)) 



inf inf /' j^\dA'{j{s))\^ds 

/ii,/i2G{0}x5O(2) -f = (x{.),Rn(-)) e ((0,1), SE (3)), ViGl3 4 5l 

7(0) = e^2>7(l) = S^Sl^l, V ^ ' ' -r 



(167) 



where we set D'^'^ = P ^, s spatial arc-length, and D^^ = D^^ = 1, so that the equivalent problem on R'^ (recall 
Eq. (pl6|)) is: 



6Zh-((y,n),(/,nO)= inf 

x(-) e c°°((o,l),r3)^ 

x(0) = 0, x(0) = e^, 
xiL) = R^iy' -y), 



(168) 



where s, L > 0, and k,{s) are respectively spatial arclength, total length, and curvature of the spatial part of the curve. 



Theorem 6 The stationary curves s x{s) G that minimize 

L 

J 7/^2^^ ds 



with arclength parameter s G [0, L], curvature k, and (3 > and free length L > 0, subject to boundary conditions 
x{0) = Xq = e M^, x(0) = no = ez ^ S'^ andx{L) = xi G M ^, x{L) = m e S'^ coincide with the stationary 
curves s {x{s)^x{s)) G x 5*^ of the variational problem (116) on M.^ x S'^ and they satisfy 



d_ 
dt 



6 

/ = with = VTRPT^ adt + V KO' 



with ds = adt, a = for all horizontal curve t Nt perturbations on the extended 15 -dimensional manifold 

Z = SE{3) X (R+)3 X {TlsE{3))y. 
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Along the the minimizing curves 7 = />: , z^: , a, Ai, . . . , Ae) in Z the Lagrange -multipliers are given by 

Ai(s) = -A5(s), A2(s) = A4(s), A3(s) = 13^1 - {X^isW - iH^W, 
Xi{s) = z^s), A5(s) = z'^is), Ae = 0, 

expressed in the normalized curvature 

z{s) = {zHs), z\s)) := -=^==(^i(.), ^\s)) = cosh(/3.) z(0) + ^^^^ z^(0), 

with K — Y^^(^^P^r7^^. They satisfy the following preservation-laws 

(Ai)^ + (A^)^ + (Aa)^ = , (co-adjoint orbits) , 
r'(A3)' + (A4)' + (A5)^ = l, 

with c > constant. The SE{3)-part s g{s) of the stationary curves in Z is obtained by means of 

A(s) m{g-'){s) = A(0) m{g-^){0) = X{L) m{g-'){L), (169) 

with row-vector A = (Ai, . . . , Ae) and where m : SE{3) ^ R^^^ is the matrix group representation given by 

ni{g) = p '^'r)' for all g = {x, R) e SE{3), (170) 

with ax G R^^^ such that axy = x x y. The curvature and torsion magnitude of (the spatial part of) these stationary 
curves are given by 

n{s) = andris) = 



with z{s) - ^/ + = ||z(s)|| with constant zo ^ Zq G R^, L > such that 

x{L) = Xi andx{L) = rii. 

The spatial part s x{s) of the curves s g{s) can be obtained by integrating the Frenet formulas 

tis) \ / k{s) \ / r(s) 
Ms) = t{s) N{s) 

B{s) J V -^(^) / V 

where again 5 G [0, L] denotes arclength parametrization of the spatial part of the curve with initial condition 

m \ f , 

MO) = ^ ((^o)2^x - {zo)iey) 
^(0) / \ ^((^0)1^0^ + ^0)2^2;) 

If these curves are a global minimizer then they coincide with the geodesies on the sub-Riemannian manifold 
(6'£^(3), d^^, d^^, d^^) subject to left-invariant metric 

= P^dA^ d^^ + dA^ dA"^ + d^^ d^^ 

Proof We follow the general theory on optimizing Lagrangians on contact manifolds, (Ml [131 and we set the 1-form 
onZ 

6 
i=l 
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where we must first find the Lagrange-multipHers. Now suppose we have a 1 -parameter family of Legendre sub- 
manifolds {A^tjtGM within Z, this corresponds to a one parameter family of horizontal vector fields on {SE{3)). Then 
compute the variation of the integrated Lagrangian-form along Nt : 

s i * = i = L, > + /„. = 1 

where we used the well-known Stokes Theorem d(^J?/^) = ^^^^ = and the formula for Lie derivatives 
of volume forms along vector fields jCxA = X\dA + d{X\A). Consequently, the optimal/characteristic curves are 
entirely determined by 

y(5)J di^l^is) =0 for all 5 >0. 
But this by definition of the insert operator J means that di/j]^ (j^v) = for all v G T{Z), or equivalently formulated 

(^J dV^I^)(7 ) = for all v e T{Z). (171) 
For further details see |[T3ll . So in order to find the optimal trajectories one has to integrate the Pfaffian system 

v\dilj = Oforall V G T(Z), (172) 

consisting of 15 vanishing dual forms. 
Now by Cartan's structural formula we have 

6 

d(d^^) = - ^ 4d^^ dA^, 



where cf^ are the structure constants of the Lie-algebra C{SE(^)), recall Eq. 28) and using this crucial identity one 
finds after some straightforward (but intense) computation the following explicitrorm of \112)\ 







= 0, 




dxMi> = 


dA^ 


= 0, 






dyl3 


- adt = 0, 




dxAd^ = 


d^4 


— K^ffdt = 


0, 


9AjdV = 




— K'^adt = 


0, 


dxM'^ = 


dA^ 


= 0, 





a^jd^ = (v4RFTF - A3 - - X5K^)dt = 0, 
a„ijdv = («Vv /||kP + /?2 _ x^)adt = 0, 
d,2\d^ = (/tVvf^pTF - A5)(7dt = 0, 

-^1 J d^ = dAi - A2d^6 + ^gd^S ..3, - - 

-A2\d^p = dX2 + Aid^e _ ^ 0^ ^^j^^ ^ ^^4 _ ^1^(1^ ^ 0^ U /-^^ 

-ylgjd^/' = dAs - XidA^ + X2dA'^ = 0, 
-Ai\d^p = dXi - X^dA^ + XedA^ - Aad^^ _^ ^^,^^2 ^ g, 
-A^jd^ = dAs - Agd^-* + A4d^6 _ ^^dA^ + XidA^ = 0, 
-Aejd^ = dAe - A4d^5 + Agd^^ _ ^^dA^ + Aad^^ = 0, 

where in the last three Pfaffian forms in the left column one may set = dA^ = 0. Now by Noether's Theorem 

the characteristic curves in Z are contained within the co-adjoint orbits, cf. L14J . As the adjoint action is given by 
push-forward of conjugation Ad((;) = (LgRg-i )^, 

3 3 

i=l i=l 

with (c^)^ + (c^)^ + (c^)^ = 1, where we used Ai\g = {Lg)^Ai and as the coadjoint action given by 

6 6 

{CoAd{g)^X,dA\^c^Aj) = {^X,dA\Ad{g-')^c^Aj) 

i=l j 1=1 j 

acts transitively on the span by the dual angular generators we see that SE{3) co-adjoint orbits are given by 

(Xif + {X2f + (Xsf = 



69 



C > 0. Indeed it directly follows from ( 172) that J2 ^jdA^ = 0. The first three Pfaffian forms in ( 172 ( yield 



A3 
A4 

As 



= , A2 — A4, 
Ai = — A5 



(174) 



so that another preservation law is given by 

r'(A3)' + (A4)' + (A5)' = l, 



22) that 



and as \q = X^hz'^ — X^hz^ = 0, Ae is constant. Now since = is an extra constraint (recall Figure 
allows us to choose convenient representatives in the cosets, the minimizers on the quotient space do not depend on 
this constraint. Therefor^ we can set Ae = and we find 



2.^2 



As + Ai = ^ As + Ai = As 



As - /^^As = 0, 



so that As = /3^As. Analogously one finds A4 = /3^A4, so similar to the {SE{2)^ d^^)-case one has 



z{s) = /3^z(s) , with z{s) 



1 



(175) 



for the normalized curvature z(s). Furthermore we find 

Ai = —XsK:^ and A2 = Xstz^. 

So again in contrast to the common and well-known elastica the geodesies do not involve special functions and we 
have 

sinh(/35) , 



z{s) = cosh(/35) z(0) + 



/3 



(176) 



Next we compute curvature magnitude and torsion magnitude. Curvature magnitude is given by 



In order to compute the torsion, we compute the normal and binormal from the tangent 

T(s) = AU- 

When taking derivatives one should use covariant derivatives, recall Theorem [3] to compute 

d 



ds 



More explicitly this yields 



£ ^3|^(,) = ^l|^(,) - >^Hs) , 

1^ AUf,) = ^2|^(,)-K2(s) 



(177) 



7(s) 



Comparison to an alternative derivation generalizing the approach in f54l on to confirmed this choice. 
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along the horizontal optimal curves in SE{3), where we used ( 165 ). Using ( 177 ) we find 



N{s) = 
N^(.) = 



^ ^ Ai 



{k'^{s)Ai - K^A2) 
Ai 



his 



A2\ 



7(5) 



B(.) 



-;{k^{s)Ai^ k'^A2) 



Now the Frenet- frame satisfies N = — /t:T + rB so curvature magnitude and torsion magnitude are given by 



n = V(/^i)2 + (a^2)2 and r 

So if we write ||z|| = we see that 

ii^ii 



^1^2 _ ^2^1 



2\2 • 



and r 



det(zo I t'q) 



(178) 



(179) 



where zq = z(0) G and Zq = z'(0) G and where we note that and z'^ satisfy the same linear ODE with 
(constant) Wronskian 

z''{s)z^{s) - z''{s)z^{s) = ^'(O)i'(O) - ^'(O)i^(O) = det(zo | 1!^). 

We verified our results by extending Mumford's approach to elastica on to curves minimizing ^/k'^{s) + 
with L free on R^. The details will follow in future work. For now we mention that such approach finally results in 
the following non-linear ODE's for n and r: 



•r2)||z||and 2tM\' ^ r' 



0. 



(180) 



It can be verified that the linear ODE's Eq. \\15\ and Eq. (179) indeed imply Eq. ( 180 ). So, besides of all the 
preservation-laws, we see another advantage of the symplectic differential geometrical approach on the extended 
manifold 15 dimensional manifold Z over the more basic approach on R^. Another advantage of the symplectic 
differential geometrical approach is the Frenet frame integration. Basically, one has to integrate 





n{s) 
-k{s) t{s) 
-t{s) 

where again 5 G [0, L] denotes arclength parametrization of the spatial part of the curve with initial condition 



^ ((^0)26:, - {Z0)l^y) 

^ ((^o)ie:, + (^0)26^;) 



(181) 




where we recall that the angles in the second chart (/3o, 70) can be computed from no by means of Eq. ( 111 ). Now 
everything is expressed in the five constants 

L,zo = ((2;o)i,(^o)2)^andzo = (4)2)^, 

which yield five degrees of freedom needed to meet the boundary condition 

x(L) = xi G R^ and x = ni G 5^ 
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Now in order to investigate on how the length L and the integration constants zq and Zq depend on the boundary con- 
dition, we follow a similar approach as in (IST, App.A] where we managed to obtain full control over the equivalenj^ 
problem in SE(2). To this end we re-express the last 6 Pfaffian forms in ( [173| ) as 



dA 



-\m{g) ^dm{g), 



(182) 



where we recall ( 170), where m : SE{3) 



3)6x6 



forms a (uncommon) group representation as we will show next. 



Taking the outer product y cr^y = x x y is tensorial with respect to rotations and ctr^ = Rcr^R ^ indeed implies 

m{gig2) = m{gi)m{g2) and m(e) / 



for all ^1, ^2 ^ SE{3). The identity ( 182| ) follows from Eq.( |173| ) by means of 

m{g)~^dm{g) = 



R-^dR R-^da^R 
R-^dR 



R ^dR cr^-idx 

R-^dR 



and recall from Eq. (33 ) that R ^dx = (d^^, d^^, d^^)^ and where the matrix g ^dg contains 
recall Theorem [3| as elements. Now from ( |182| ) we deduce that 

d{Xm{g)-') = Xd{m{g)-') + {dX){m{g))-' = {X{m{g))-' dm{g) + dA) {m{g))-' = 

and the result dl69|) follows. 



□ 



Theorem 7 The Frenet system \181) for the stationary curves can be integrated analytically. By means ofEq. \169\ 
we find 

m{g) = hQ^m{g), g = (x, y, z, R), (183) 

where h^^ e SE{3) is given by 



— T — T 

rI 



{m{xo,Ro)) ^=m{{xo,Ro) ^) 



withxo = (xo, Vo^zq) e and Rq G 5*0 (3) given by 



xo = (free), 

^0 — r.mWLW ' 



(184) 



(185) 



with Wronskian W = Zi{s)z2{s) — Z2{s)zi{s) = ^i(0)i2(0) — ^2(0)^1(0) <^nd 



Ro — 



cf3 



-(4)2 

117II1 



112 



K\\ 114 II 





114 II 



withzQ = {{zq)i, (4)2) if^o 7^ ^- case Zq = we have 



Xq 



1 





(^o)i 
(^0)2 



andRo ^e^,f • 



Furthermore, we have 



g{s) = {i{s),R{s)) = {x{s),y{s),~z{s),R{s)) 



(186) 



(187) 



^^In 1 31 App.A] we had to choose zq and z'q and L such that x(L) = xi G and x(L) = ni ^ . 
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for all s G [0, L] with 



x{s) = m -l^^[E [{Ps + f )z, m)-E ((f )z, m)) , 

S A{s')ds' { n 11^(^)11^ 

{y{s), z{s))^ = eo (^(0), ^(0))^ = Re { {zq - i^o)e zM.i(.)-^wVl-o-^^-^m 



where the elliptic integral of the second kind is given by E{(j)^ m) = j \/l — m sin^ OdO with m = and with 



^ ^ \\zo-r'zMzo+r'z^o\\ < ^^Y/z (/) = log and with (commuting) 



A{s) 



z{s)-z{s) 



z{s)-z{s) 



and with R{s) in S0{3) such that 



R{s)e,=n{s) = {S:{s),{A{s){y{s),z{s)ff) 



Proof With respect to Eq. ( |183| l,( |184[ ), ( |185| i, ( |186| l and ( |187[ > we note that by means of Eq. ( |169| l we have 

A(m(5))-i = /X 4^ A = /xm(5) = ^ih^^m{~g) = \{Q){m{g{m~^ ^{g) = \{Q)h^^m{g) 

as 5(0) = (0, 1) G SE{3) implies that its corresponding matrix representation equals the 6 x 6-identity matrix 
m{g{0)) = Iq. Then we choose h^^ such that 

Ai/io' =A(0)/ioi = (c/?,0,0,- — ,0,0) 



with /i^ = A^(0), i = 1, . . . , 6, which follows from the previous theorem. As a result Hq now follows from 

Rq 



1 ) 






C/3 























w 

c/3 




/i4 



















V y 




V 





J 



cp\ _ / (zo)i 
xo X I j +i^o I (^0)2 








ci3 






(188) 



Next we shall derive the ODE for z). We parameterize by R = R^^^^R^ ^Re^,a ^ S0{3) and substitution of this 
parametrization in 



W 



yields the system 



and the system 



Ai 

i?| A2 

A3 



/XfTAf A3 
c/3 

-A3 v^AfTAf 



c/3 





- ^2 





A2 







1 / 



A4 

As 




(189) 



_w 

cl3 

-c/3z 
c/Sy 



(190) 
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where we used P = arg(^Af + -\- iXs) and a = arg(Ai —1X2) that follow by Eq. (189). Taking the derivative 
with respect to spatial arc-length s in system ( |19Q| ) yields 

^R{X^^Y^ = R{X^^Y ^ R{X^^Y = -I y withA^2) = (A4,A5,0), 



(is 

where we stress that 

R^dR = cr((i^4^d^5^o) = ^(/■cids,/^2ds,0) =^ 

^(A(2)f = ^a(,i,,2,o)T 1^ A5 j = ^ 
As a result, since (Ai)^ + (-^2)^ = ||z||, we obtain 



y 



^^^and^=^^^, (191) 



whereas ( |19Q| ) yields 



from which we deduce that 



1^ COS7 ^^g^ sm7 = -cf3z, 



. ^ , H^V'c2^2_||^||2 ^ 

sm7 + ^,.||,|| COS7 = cpy 



(Ms))-' ( ] = ( ^S;\ ] ^ ( %i ] = A{s) ( y^'^ 



z{s) J \ z{s) J \ z{s) J \ ^(^) 

As the matrices A{s) share the same eigen-vectors ((z, 1)^ and (— 1)^) they commute and we have 

( f{s) )=e!SMr)^r ^ m ^ ^ 1.(^(0) +,,~(Q))e/o %(.)...)i^w^!:(.)/Ce.) Ci. ^ j 

where the latter two vectors are each other's conjugate. Finally, x{s) follows by integration 
x{s) =x(0) + /^dr 



x(0) + i/Vl-||z(r)Pdr 





x(0) + / a/1 - ci cosh(2/3r) - C2 sinh(2/3r) dr 

cV2 Q 



withci = -^Y^^dlzolP - ||zo|P), C2 = j^^^^^q ■zo,m= which can be expressed as 

c^/2 







: i(0) - + f )i, m) - i? ((f )i,m)) 



where ^ = ^ log ^^^^^ , v = i^r, j = y/ (ci)^ — (02)^ and where we note that 



2 ^^fc> C1-C2 

1 11 . Ml. . . 11 /,, ||2 I II I p-1- \\2\ 1 — 



zo-/3-izo||||zo + /3-izo|| < ^^-(||zo-/3-izo|p + ||zo + /3-izo|p) = <l.n 



1+C2" " ' " ' ""-1+C22^" ' " " ' 1+C 
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Figure 23: Example of a geodesic in x S'^. Left: s x(s). Right: s ^ x(5). Parameter settings = jq, 
zo = (^,0)^, z[) = -/3zo, (so c = 1 and 1^ = 0), length L = 30, x(0) = (0,0,1)^, x(0) = (0,0,0), the 
rotation and translation needed to map the curve x onto the geodesic x is gi ven b y h^^ = (m(xo, i?o))~^- Note that 



J" _ _ 

x{s) = Rq {x{s) — xo), with Xo = (0, 0, —5) and Rq G 50(3) given by Eq. ( 186). In this example both the curvature 



K, and z are always aligned with the y-axis which explains that the curve stays within the yz-pla.m. 

Corollary 3 The solution curves stay in a plane iffzo and Zq are linear dependent, i.e. W = det(zo|zo) = 0. In the 
special case Zq = —Pzo, (c = 1 and W = 0) we find 

x{s) = m + « + i (^/r^w - v^r^w^ + log (i±^^^^g^ 

In the case where solution curves stay in the plane, i.e. W = {), with general c, we find 



i(s)=i(0) + i/Vl-||z(r)PdT, 
y{s) = 



and 

~z{s) = 



sgn{so-s)~z{0)^S\ f^' ^ > ^ 



where z{s) is given by Eq. {176) and = ^(/n(||zo|p - /3^||zo|P) - lri{\\zo + /^ZolP)) so that \\z{so)\\ = 0. These 
solutions indeed coincide with the geodesies on the contact manifold {SE{2) = x S*^, — sin 6dx + cos Ody) where 
elements are parameterized by (x, e*^), cf. [34 App.A]. 



See Figure 23 for an example of a geodesic with W = and c = 1. 
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